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i t rat, if ication,  i.e.,  density  variations  in  a  reservoir,  occurs  due  to 
temperature  variations  as  a  result  of  surface  heat  exchange  and  plays  an  impor¬ 
tant  role  in  determining  the  water  quality  of  a  reservoir.  This  role  is  deter¬ 
mined  through  the  influence  of  density  variations  on  the  movement  of  w  iter  in 
the  reservoir.  Therefore,  the  primary  objective  of  a  prediction  of  stratified 
flow  hydrodynamics  in  reservoirs  is  to  enable  scientists  to  compute  temperatures 
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distributions  and  water  transports  insofar  as  they  affect  various  water  quality 
parameters. 

One  objective  of  the  Environmental  &  Water  Quality  Operational  Study 
(EWQOS)  program  of  the  U.  S.  Army  Corps  of  fcngineers  is  to  provide  District 
and  Division  offices  with  a  tool  for  predicting  reservoir  hydrodynamics  over 
periods  of  time  extending  from  the  initial  sptup  of  thermal  stratification  in 
the  spring  through  its  breakup  in  the  fall.  VSuch  a  predictive  technique  will 
subsequently  be  used  in  the  prediction  of  wat|pr  quality  parameters.  /  An  impor¬ 
tant  tool  that  provides  a  relatively  low  cost  highly  flexible  model  »of  the 
hydrodynamics  of  a  water  body  is  a  numerical  model. 

Under  an  EWQOS* work  uniX^T^pth  two-  and  three-dimensional,  unsteady, 
variable  density,  heat-conducting'models  have  been  investigated  during  the 
past  year.  This  investigation  has  centered  around  an  analysis  of  both  the 
mathematical  and  numerical  bases  of  individual  models  as  well  as  their  ability 
to  simulate  a  density  underflow^n  the  Generalized  Reservoir  Hydrodynamics  (GRH) 
flume  located  at  the  U.  S.  Army  ^Engineer  Waterways  Experiment  Station  (WES). 

A  discussion  of  the  limitations  land  relative  advantages  of  the  various  models 
is  presented  along  with  results  tf  the  GRH  flume  applications. 

The  general  conclusion  is  that  a  two-dimensional  laterally  averaged  model 
developed  by  Edinger  and  Buchak  offers  the  most  promise  of  providing  the  Corps 
with  a  computationally  efficient  and  accurate  multidimensional  reservoir  hydro- 
dynamic  model.  It  does  not  appear  that  any  three-dimensional  models  that  allow 
for  economical  long-term  simulations  have  been  developed. 
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CONVERSION  FACTORS,  U.  S.  CUSTOMARY  TO  METRIC  (Si) 
UNITS  OF  MEASUREMENT 


U.  S.  customary  units  of  measurement  can  be  converted  to  metric  (SI) 
units  as  follows: 


_ Multiply _ 

cubic  feet  per  second 
Fahrenheit  degrees 
feet 

feet  per  second 

pounds  (mass)  per  cubic  foot 

square  feet 

square  feet  per  second 


-  by _ 

0.02831685 

5/9 

O.30U8 

O.30U8 

0.016018U6 

0.09290301* 

0.0929030U 


_ To  Obtain _ 

cubic  metres  per  second 
Celsius  degrees  or  Kelvins* 
metres 

metres  per  second 

grams  per  cubic  centimetre 

square  metres 

square  metres  per  second 


*  To  obtain  Celsius  (C)  temperature  readings  from  Fahrenheit  (F)  read¬ 
ings,  use  the  following  formula:  C  =  ( 5/9) (F  -  32) .  To  obtain 
Kelvin  (K)  readings,  use:  K  =  (5/9)(F  -  32)  +  273.15. 
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A  REVIEW  OF  NUMERICAL  RESERVOIR 
HYDRODYNAMIC  MODELING 


PART  I :  INTRODUCTION 


Reservoir  Straiification  and  Its  Importance 

1.  As  the  population  of  the  United  States  has  increased  over  the 
past  few  decades,  there  has  been  a  corresponding  increase  in  the  demand 
on  water  supplies.  To  help  meet  this  demand,  numerous  impounding  reser¬ 
voirs  have  been  constructed.  The  impoundment  or  damming  of  a  flowing 
stream  can  significantly  affect  the  quality  of  the  water.  This  can  hap¬ 
pen  as  a  result  of  the  direct  increase  in  travel  time  required  for  water 
to  traverse  the  distance  from  the  headwater  of  the  stream  to  the  dam  as 
well  as  the  effect  that  stratification  plays  in  determining  the  quality 
of  the  water  released  from  the  reservoir.  The  relationship  between 
density  variations  and  quality  parameters  in  the  reservoir  is  a  direct 
result  of  the  influence  of  stratification  on  the  movement  and  mixing  of 
water. 

2.  Stratification  or  density  variations  in  a  reservoir  can  occur 
as  a  result  of  solute  concentrations,  suspended  solids  concentrations, 
or  temperature  variations  as  a  result  of  surface  heat  exchange.  Surface 
heat  exchange  is  a  function  of  both  short-  and  longwave  radiation  as 
well  as  surface  conduction,  evaporation,  and  precipitation.  In  the 
remainder  of  this  report,  the  term  "sti atif ication"  will  refer  to  density 
variations  due  to  thermal  effects. 

3.  At  the  beginning  of  spring,  a  reservoir  is  essentially  homoge¬ 
neous.  However,  as  the  weather  warms,  the  water  near  the  surface  also 
warms  due  to  an  exchange  of  heat  from  the  atmosphere  to  the  water  sur¬ 
face.  The  warmer  water  near  the  surface  is  then  mixed  downward,  primarily 
by  wind  action  and  diurnal  cooling.  By  late  summer,  the  reservoir  will 
attain  maximum  stratification  (Figure  l).  A  warm  upper  layer  (epilimnion) 
of  water  resides  over  the  cold  deeper  layer  (hypolimnion)  with  a  zone 
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Figure  1.  Regions  associated  with  thermal 
stratification 

between  the  two  (metalimnion)  in  which  a  large  density  gradient  exists. 
As  the  weather  cools  during  the  fall,  the  surface  temperature  decreases, 
resulting  in  denser  water  at  the  surface  and  a  correspond! np  convective 
overturning.  This  mixing  eventually  results  in  an  i sotherr.al  water  body 
that  remains  isothermal  through  the  winter,  except  during  periods  of  ice 
cover.  Such  a  cyclic  variation  of  temperature  is  demonstrated  by  the 
seasonal  temperature  profiles  presented  in  Figure  ? . 

Density  Currents 


.  The  variation  of  the  fluid  density  in  a  stratified  r«o*-rv.  i  r 
gives  rise  to  what  are  known  as  internal  density  current  s  or  s’ rat  r 
flow.  Such  flow  refers  to  motions  involving  fluid  masses  of  *  h>  ■•uve 
phase.  A  heavier  liquid  flowing  beneath  a  lighter  liquid  r  m  h>avi<r 


JUNE  9 


JULY  26 


MARCH  1 


TEMPERATURE,  °F 

Figure  2.  Temperature  profiles  of  the  west 
basin  of  Horn  Lake,  B.  C.,  during  I960 
(taken  from  Slotta  et  al.  1969) 


gas  moving  under  a  lighter  gas  will  be  subject  to  gravitational  effects 
that  depend  upon  the  differences  between  the  two  densities.  Such  flows 
can  be  extremely  important.  Slotta  et  al.  (1969)  discusses  an  example 
of  a  density  current  in  the  Watts  Bar  reservoir  of  the  Tennessee  Valley 
Authority  (TVA)  system  in  which  a  coldwater  density  flow  moves  over 
13  miles*  upstream  into  the  warmer  waters  of  one  of  the  arms  of  the 


*  A  table  of  factors  for  converting  U.  S.  customary  units  of  measure 
ment  to  metric  (Si)  units  is  presented  on  page  U. 
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reservoir.  This  bottom  density  current  flows  past  a  sever  outlet  as 
well  as  the  outfall  from  a  large  paper  mill.  At  or.e  time,  sewage  and 
mill  waste  were  discharged  into  the  coldwater  current  and  carried  up¬ 
stream  to  the  intake  of  a  water  plant.  The  situation  since  teen 
corrected  by  use  of  a  variable  level  outfall  for  the  sewage  and  mill 
waste.  The  internal  motions  in  reservoirs  due  to  temperature  variations 
or  perhaps  due  to  the  inflow  of  sediment-laden  streams  plus  the  under¬ 
standing  and  control  of  salinity  intrusion  in  tidal,  estuaries  are  amen r 
the  most  challenging  of  present-day  problems  dealing  witn  strati  fie  d 
flow. 


Relationship  of  Water  Quality  to  Hydrodynamics 

5.  The  primary  objective  of  a  prediction  cf  stratified  flew  in 
reservoirs  is  to  enable  scientists  to  compute  temperature  distributions 
and  water  transports  insofar  as  they  affect  various  water  quality 
parameters.  While  the  process  of  heat  transfer  in  bodies  cf  water  is 
nothing  new,  the  prediction  cf  the  resulting  multidimensional  flow 
phenomena  in  a  reservoir  for  varying  stream  inputs  as  well  as  dam  dis¬ 
charges  from  varying  levels  is  extremely  difficult. 

6.  A  substance  (either  chemical  or  biological)  disperses  through 
a  reservoir  by  convection  and  turbulent  diffusion.  In  addition,  the 
substance  is  also  acted  upon  by  various  chemical,  biological,  and  physi¬ 
cal  processes.  An  understanding  of  both  the  dispersion  and  the  chemical 
and  biological  processes  is  essential  in  any  prediction  cf  water  quality, 
which  is  the  ultimate  goal  to  be  sought,  although,  not  the  coal  of  this 
study.  It  should  be  clear  then  that  a  problem  of  such  scope  calls  for 

a  cooperative  effort  cf  a  wide  variety  of  scientific  disciplines  rang¬ 
ing  from  meteorology,  hydrology,  and  hydrodynami cs  to  chemistry  and 
biology. 


Hydrodynamic  Predictive  Techniques 

7.  In  an  attempt  to  predict  the  hydrodynamics  of  a  reservoir,  one 
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or  perhaps  a  combination  of  three  approaches  may  be  taken — field 
investigations,  physical  models,  and  mathematical  or  numerical  models. 
Field  investigations  may  reveal  what  presently  occurs  in  a  water  sys¬ 
tem,  but  cannot  predict  what  will  result  from  changes  due  to  new  inputs 
to  the  system.  In  addition,  field  investigations  are  usually  relatively 
expensive.  Depending  upon  their  complexity,  e.g.,  large  models  of 
river  basins,  estuaries,  and  bays,  physical  models  can  require  signifi¬ 
cant  investments  of  capital,  long  construction  times,  and  long  test 
periods.  However,  physical  models  of  reservoirs  to  address  problems 
such  as  near-field  inflow  selective  withdrawal  and  pumpback  character¬ 
istics  of  specific  outlet  structures  and  geometries  are  far  less  expen¬ 
sive  to  construct  and  operate.  Depending  upon  approximations  made  to 
the  governing  equations  of  motion  and  the  solution  technique  employed 
to  solve  the  equations,  numerical  models  can  often  provide  relatively 
low  cost  and  highly  flexible  models.  However,  it  should  be  remembered 
that,  as  with  many  physical  models,  numerical  models  must  be  calibrated 
and  verified  before  confidence  can  be  placed  in  results  obtained  from 
them.  Data  from  both  field  studies  and  physical  model  tests  are  used  to 
assess  the  reasonableness  of  numerical  predictions  and  to  aid  in  the 
further  development  of  mathematical  descriptions  of  observed  physical 
processes.  However,  the  steady  advances  in  computer  technology  over 
the  past  two  decades  (Table  l)  indicate  the  potential  for  even  greater 
economical  use  of  more  widely  applicable  numerical  models  in  the  future. 
In  a  practical  sense,  a  combination  of  the  most  desirable  features  of 
both  physical  and  numerical  models  will  probably  continue  to  provide  the 
best  approach  for  solving  most  hydrodynamic  problems. 


Types  of  Numerical  Models 

8.  Numerical  hydrodynamic  models  can  differ  widely,  depending 
upon  such  things  as  the  solution  technique  applied  to  the  governing 
differential  equations  representing  the  physical  processes,  the  assump¬ 
tions  made  in  the  derivation  of  the  governing  equations,  whether  the 
phenomena  are  steady  or  time-varying,  and  the  spatial  dimensionality 


considered,  with  perhaps  the  spatial  dimensionality  being  the  most 
commonly  used  delineator. 

9.  One-,  two-,  and  fully  three-dimensional  numerical  models  that 
provide  for  the  simultaneous  solution  of  the  coupled  turbulent  velocity 
field  and  the  temperature  field,  subject  to  varying  boundary  conditions, 
exist  and  are  applicable  in  varying  degrees  to  the  problem  of  predicting 
stratified  reservoir  liydrodynamics  for  use  in  developing  water  quality 
predictions.  Perhaps  the  earliest  work  in  which  computations  for  the 
fluid  density  and  the  flow  were  coupled  was  the  work  of  Welch  et  al . 
(1966)  in  tiie  development  of  the  two-dimensional  Marker  and  Cell  code 
commonly  referred  to  as  MAC.  Paralleling  the  development  of  MAC  and 
the  MAC-related  codes,  e.g.  ,  Chan  and  Street's  (1970)  SUMMAC,  Slotta 
et  al.'s  (1969)  NUMAC,  etc.,  have  been  a  host  of  models  that  might  be 
described  as  control  volume  models.  With  this  method,  a  reservoir  is 
divided  into  a  number  of  horizontal  layers  extending  over  its  breadth 
and  length.  Homogeneity  is  then  assumed  in  each  layer.  The  result  is 
a  one -dimensional  model  with  variations  allowed  only  in  the  vertical. 
Governing  differential  equations  are  obtained  by  applying  mass,  momentum, 
and  heat  balances  for  the  control  layers.  Inflow  and  outflow  boundaries 
can  be  included  quite  easily  in  such  models.  Parker  et  al.  (1975)  re¬ 
viewed  such  one-dimensional  reservoir  models  and  concluded  that  such 
models  could  be  applied  to  larger,  deep  reservoirs  where  horizontal  flow 
has  minimal  impact  on  the  vertical  density  structure.  The  primary  ad¬ 
vantage  of  such  a  model  is  its  ability  to  resolve  long-term  or  seasonal 
temperature  profiles  economically.  However,  it  must  be  noted  that  such 
one-dimensional  models  are  not  applicable  to  the  problem  with  which  this 
study  is  concerned — predicting  multidimensional  flow  fields  within 
stratified  reservoirs  for  quality  predictions. 

10.  Both  two-  and  three-dimensional  hydrodynamic  models  are  dis¬ 
cussed  in  detail  in  succeeding  sections.  Some  of  the  models  investi¬ 
gated,  such  as  the  two-dimensional  models  of  Edinger  and  Buchak  (1979), 
Waldrop  and  Farmer  (1976),  and  Roberts  and  Street  (1975),  are  directly 
applicable  to  reservoirs,  although  in  varying  degrees;  while  others, 
such  as  the  two-dimensional  depth-averaged  models  of  Leendertse  (1967), 
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Masch  et  al.  (1969),  and  Reid  and  Bodine  (1968),  have  no  applicability 
to  the  modeling  of  internal  flows  in  stratified  reservoirs  other  than 
perhaps  in  the  numerical  techniques  employed. 

11.  Three-dimensional  hydrodynamic  models  have  only  recently 
been  developed  to  the  state  where  application  to  complex  geometries  with 
reasonable  resolution  for  short  simulation  periods  appears  possible; 
however,  the  cost  is  still  prohibitive  for  simulations  over  long  periods 
of  time.  Thus,  it  appears  that  if  one  is  only  interested  in  the  steady- 
state  flow  and  temperature  field  resulting  from  situations  such  as  a 
discharge  of  warm  water  from  a  power  plant,  three-dimensional  modeling 
appears  feasible.  However,  if  the  interest  lies  in  computing  reservoir 
hydrodynamics  over  a  stratification  cycle,  i.e.,  several  months,  new 
developments  in  solution  techniques  and  the  availability  of  larger  and 
faster  computers  must  be  realized  before  such  modeling  becomes  practical. 


Purpose  and  Scope 


12.  The  need  for  a  predictive  capability — numerical  models  in 
the  area  of  stratified  reservoir  hydrodynamics — has  been  firmly  es¬ 
tablished.  The  purpose  of  the  study  described  herein  then  is  to  select 
the  most  applicable  existing  models  and  to  provide  recommendations  for 
additional  developmental  work  needed  to  meet  that  need.  Because  of  the 
nature  of  the  problem  to  be  addressed,  selected  models  must  have  the 
capability  of  handling  free  surface  variable  density  flows  that  are 
time-dependent.  PARTS  II  and  III  present  a  detailed  discussion  of  the 
theoretical  basis  and  corresponding  numerical  techniques  that  are  common 
to  all  numerical  hydrodynamic  models.  Three-dimensional  hydrodynamic 
models  are  discussed  in  PART  IV,  and  two-dimensional  hydrodynamic  models 
are  discussed  in  PART  V.  In  addition  to  an  investigation  of  the  theo¬ 
retical  limitations  of  various  models,  a  limited  attempt  at  analyzing 
the  actual  performance  of  several  models  has  been  made.  This  was  accom¬ 
plished  through  model  applications  to  a  coldwater  underflow  in  the  U.  R. 
Army  Engineer  Waterways  Experiment  Station  (WES)  Hydraulics  Laboratory’s 
General  Reservoir  Hydrodynamics  (GRH)  flume.  These  results  are  presented 


1? 


in  PART  VI.  Finally,  conclusions  of  this  study  and  recommendations  for 
additional  developmental  work  needed  to  provide  the  Corps  of  Engineers 
with  computer  models  with  the  potential  to  provide  a  predictive  capabil¬ 
ity  in  the  area  of  reservoir  hydrodynamics  in  an  accurate  and  economical 
fashion  are  presented  in  PART  VII. 
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PART  II:  MATHEMATICAL  DISCUSSIONS 


Basic  Equations  and  Approximations 


13.  The  Navier-Stokes  equations  express  the  conservation  of  mass 
and  momentum  of  a  flow  field  and  are  the  basic  governing  equations  for 
the  solution  of  any  fluid  dynamics  problem.  These  equations  written  in 
tensor  notation  are* 


3 '  ^pu i 

Cont  inuity  :  -rp  +  — -  =  0 


3x . 
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(1) 


Momentum: 


3pu. 
_ 1_ 
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Sr  .  fi.pu,  +  — ^ 
:jk  j  k  9x 
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where 


p  =  water  density 
t.  =  t  ime 

u.  =  tensor  notation  for  velocity 

x.  =  tensor  notation  of  spatial  coordinate 

g .  =  acceleration  of  gravity 


1  '  K 


-  cyclic  tensor 
=  Coriolis  parameter 


T.,  =  laminar  stress  tensor 

1.) 

u  =  molecular  eddy  viscosity 

6.  ,  =  Kronecker  delta 
ij 


and  where 


represents  the  viscous  molecular  stress  arising  as  a  result  of  the  con¬ 
tinuum  approach.  It  will  be  recalled  from  tensor  theory  that  repeated 


*  For  convenience,  symbols  are  listed  and  defined  in  the  Notation 
(Appendix  C). 
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indices  imply  a  summation  and  also  that 
the  cyclic  tensor  defined  as 


c...  in  the  Coriolis  term  is 
ijk 


e. .  =  1  ,  for  an  even  permutation  of  ijk 

1  j  h 

=  -1  ,  for  an  odd  permutation  of  i j k 
=  0  ,  otherwise 


For  example,  =  e  =  z  =  1  ;  whereas,  t  =  e  =  =  -1 

and  the  Kronecker  delta  6. ,  is  defined  as 

ij 

6ij  *  1  •  if  i  =  J 
=  0  ,  otherwise 

In  addition  to  the  above  equations,  a  conservation  of  energy  equation 
must  also  be  written  for  fluid  flow  problems  with  thermal  effects.  With 
the  assumption  of  a  constant  specific  heat  and  with  the  neglection  of 
viscous  dissipative  effects,  one  can  write  the  energy  equation  as  the 
following  transport  equation  for  temperature  T  : 


where  D.  is  equal  to  the  diffusivity  coefficient.  This  equation 
^  J 

states  that  the  temperature  can  change  as  a  result  of  advection  by  the 
flow  field,  molecular  diffusion,  and  the  actions  of  any  sources  and 
sinks  of  heat.  As  a  matter  of  fact,  this  same  equation  applies  to  the 
transport  of  any  constituent  0  ,  where  0  would  replace  T  in  the 
equation  and  appropriate  sources  and  sinks  and  boundary  conditions 
would  be  prescribed.  For  example,  in  the  numerical  modeling  of  the 
hydrodynamics  of  an  estuary,  a  similar  transport  equation  for  the  salin¬ 
ity  would  be  required. 

l4.  One  additional  equation  remains  to  be  written  in  order  to 
close  the  system.  An  equation  of  state  expressing  the  density  as  a 
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function  of  the  temperature  and  pressure  (and  salinity  in  estuarine 
modeling)  must  be  employed: 


Equation  of  State:  p  =  p(T,P)  (h) 

With  the  closure  of  the  system,  there  exist  six  equations  to  be  solved 

for  the  six  unknowns — density  p  ;  three  velocity  components  u  ,  v  ,  | 

w  ;  pressure  P  ;  and  temperature  T  . 

Time  averaging  for  turbulent  flow 

15*  The  above  equations  written  with  molecular  values  of  viscosity  , 

and  diffusivity  are  only  applicable  in  a  practical  sense  tc  laminar 
flow  fields  where  the  flow  and  thermodynamic  variables  do  not  exhibit 

random  irregular  fluctuations  in  time.  However,  most  fluids  in  motion  ! 

exhibit  such  fluctuations  and  are  referred  to  as  "turbulent  flews."  ^ 

16.  Following  Reynolds,  the  approach  normally  taken  tc  make  the  ,  j 

equations  applicable  to  turbulent  flows  is  to  assume  that  the  dependent 
variables  are  composed  of  an  average  time-varying  component  plus  a  small  | 

randomly  varying  component  about  the  average  value.  This  is  illustrated 
below. 


Thus,  one  writes 

ui(x,y,z,t)  =  u.(x,y,z,t)  +  uj(x,y,z,t.) 
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where 


t+At/2 

u^x.y.z.t)  dt 

and 

t+At/2 

u!(x,y,z,t)  dt  =  0 

u!  =  deviation  between  instantaneous  velocity  and  time-averaged 
velocity 

=  time-averaged  velocity 
At  =  time  step 

With  all  the  dependent  variables  written  in  the  form  above,  substitution 
into  Equations  1,  2,  and  3  and  then  integration  over  the  time  increment 
At  produces  the  same  form  of  the  previous  equations,  but  now  written 
with  the  time-averaged  components  as  the  dependent  variables,  plus  the 
additional  terms 


k  / 


t-At/ 


1_ 

At 


t+At/2 

/ 

t-At/2 


u!u' 

i  i 


dt 


and 


t+At/2 

k  /  dt 

t-At/2 


where  T'  =  deviation  between  instantaneous  and  time-averaged  temperature 
IT.  The  first  term  is  referred  to  as  the  turbulent  Reynolds 
stress,  since  the  high  frequency  turbulent  fluctuations  manifest  them¬ 
selves  as  viscous  stresses  acting  on  the  average  component  of  flow. 

'Tsing  Boussinesq's  concept  of  eddy  viscosity,  the  first  term  is 
written  as 


1_ 

At 


t+At/2 

/ 

t-At/ 


'  d  W  <- 
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at 
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+ 


summation  over 
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In  analogy  with  the  laminar  flow  case,  e  is  referred  to  as  the 
turbulent  or  eddy  viscosity  tensor. 

18.  In  a  similar  fashion,  the  second  term  above,  which  arises 
from  the  time  averaging  of  the  temperature  equation,  is  commonly 
written  as 


1_ 

At 


t+At/2 

/ 

t-At/2 


T'u!  dt 

i 


where  A.  is  called  the  "eddy  diffusivity  tensor"  and  T  is  the  time- 
^  J 

averaged  temperature. 

19.  The  equations  commonly  applied  to  turbulent  flow  problems 
can  now  be  written  as 

3—  3pu. 

Continuity :  —  +  — —  =0  ( 5  ) 


3pu.  3(pu.u.) 


Momentum:  ^  ■  + 


i  J  _  _  3P 


3x, 


3P 

377  +  pgi 
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-  2t. ,,  n.pu.  +  t — 
ijk  j  k  3x 


C3u .  3u . 
377  +  3^- 


(6) 


f  +  37"  =  ^77  (Aij  fj)  +  Z  sources  -  Z  sinks  (t: 


Equation  of  State:  p  =  p(T,P) 


(8) 


t 

t 

l 


where 

p  =  time-averaged  water  density 
P  =  time-averaged  pressure 

and  where  the  assumption  has  been  made  that  the  eddy  coefficients  are 
much  larger  than  the  molecular  values;  i.e.. 
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ij 


Boussinesq  approximation 

20.  Subject  to  the  assumptions  made  in  their  derivation,  the  time- 
averaged  governing  equations  (Equations  5-8)  are  applicable  to  any  turbu¬ 
lent  fluid  dynamics  problem.  An  approximation  usually  made  when  applying 
the  equations  to  hydrodynamic  problems  was  first  pointed  out  by 
Boussinesq.  When  variations  of  temperature  are  small  (At  _<  +10°C), 
variations  in  density  will  be  less  than  one  percent.  For  example, 

Edinger  and  Buchak's  expression  relating  the  density  of  water  to  the 
water's  temperature  results  in  only  a  0.15  percent  increase  in  the  den¬ 
sity  when  decreasing  the  temperature  from  20°C  to  10°C.  These  small 
variations  can  be  ignored  in  general  with  one  exception.  The  variability 
of  density  in  the  gravitational  term  cannot  be  ignored.  Hence,  p  is 
treated  as  a  constant  in  all  places  except  the  body  force  term. 

21.  With  the  Boussinesq  approximation,  the  continuity  and  momentum 
equations  become 

3u. 

Continuity:  =0  (9) 
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Momentum:  — —  +  — = 
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3x 
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Pq  3x.  Pq  i 


2e.  n  u,  +  —  -r— 
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/ 3u^  3u 
eij\3xj  +  3xi; 


where  p  is  a  reference  water  density.  The  energy  equation  and  the 
o 

equation  of  state  are  not  affected. 

Conservative  versus  nonconservative 

22.  When  the  momentum  equations  are  written  as  Equation  10,  they 
are  known  as  the  conservative  form  of  '■-he  equations.  If  the  convective 
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term  is  expanded  and  Equation  9  is  substituted.  Equation  10  reduces  to 


3u.  _  3u.  _p  /—  v  .  _ 

3t“+Uj  3^  =  -  ^"3Z'  +  \  Jj  gi  ~  2£ijkttj\ 
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which  is  referred  to  as  the  nonconservative  form.  Analytically,  the 
two  forms  are  equivalent.  However,  in  numerical  solutions  of  flow  prob¬ 
lems,  they  are  not.  As  discussed  by  Leendertse  (1967),  a  finite  differ¬ 
ence  representation  of  Equation  11  does  not  conserve  momentum  of  the 
flow  field;  whereas,  the  identical  numerical  representation  of  Equa¬ 
tion  10  does.  As  a  result,  most  of  the  more  recent  numerical  hydro- 
dynamic  models  use  the  conservative  form  as  opposed  to  the  nonconserva¬ 
tive  form  employed  in  many  of  the  earlier  models. 

23.  An  interesting  point  is  that  in  the  laminar  form  of  the 
momentum  equations,  i.e.,  Equation  2,  when  the  assumption  of  incompres¬ 
sibility  is  made,  researchers  have  historically  neglected  that  portion 
of  the  viscous  term  that  contains  the  condition  of  incompressibility, 

2  3ui 

i.e.,  —  p  6  ,  even  though  they  may  have  retained  the  conservative 


form  of  the  convective  terms.  In  the  turbulent  form  of  the  equations, 
there  is  no  such  inconsistency,  since  all  molecular  viscous  terms  are 
neglected  due  to  the  assumption  that  the  eddy  viscosity  is  much  larger 
than  the  molecular  viscosity.  Therefore,  the  condition  of  incompressi- 

term  in  the  turbulent 


2  ?Ui 

bility  is  not  invoked  in  dropping  the  —  p  6,.  . 


j 


form  of  the  equations. 

Convective  versus  quasi-static 

2b.  An  assumption  that  is  usually  made  in  hydrodynamic  models  is 
that  vertical  accelerations  are  negligible  when  compared  to  the  gravita¬ 
tional  acceleration.  Neglecting  viscous  terms  also,  the  vertical  momen¬ 
tum  equation  reduces  to 
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(12) 


3P 


3z 


Pg 


where  z  is  the  vertical  spatial  coordinate  and  is  positive  downward. 
Equation  12,  of  course,  states  that  the  pressure  is  hydrostatic. 

25.  When  considering  the  coupling  of  the  thermodynamics  and  the 
hydrodynamics  of  a  water  body,  a  distinction  must  be  made  between  con¬ 
vective  and,  as  labeled  by  Simons  (1973),  "quasi-static”  models.  Con¬ 
vective  models  retain  the  complete  vertical  momentum  equation  and  can 
simulate  in  full  detail  the  convective  overturning  of  unstable  water 
masses,  such  as  the  upwelling  of  cells  of  warm  water  or  perhaps  the 
plunging  of  a  ccldwater  inflow.  Quasi-static  models  where  the  pressure 
is  hydrostatic  eliminate  vertical  accelerations  due  to  buoyancy  effects, 
which  precludes  the  explicit  treatment  of  free  convection  associated 
with  unstable  stratification.  Convective  overturnings  can  only  be 
handled  as  mixing  along  the  vertical  axis. 

26.  A  commonly  used  technique  is  that  of  invoking  a  large  verti¬ 
cal  diffusion  of  heat  to  counteract  such  instabilities.  This  results  in 
the  removal  of  any  unstable  stratification  the  moment  it  occurs.  Such 

a  technique  is  implemented  by  checking  the  vertical  temperature  profile 
at  each  horizontal  location  after  each  computation.  If  at  any  point 
lighter  water  lies  below  denser  water,  the  profile  is  adjusted  without 
affecting  the  total  heat  content  of  the  column. 

27.  As  will  be  discussed  in  more  detail  in  connection  with  the 
numerical  solution  of  the  governing  equations,  models  with  the  hydro¬ 
static  assumption  require  far  less  computer  time  than  the  fully  convec¬ 
tive  models. 

Spatial  averaging 

28.  A  solution  of  Equations  7,  8,  9,  and  either  10  or  11  consti¬ 
tutes  a  fully  time-varying,  three-dimensional  model  with  the  only  assump¬ 
tion  being  the  Boussinesq  approximation.  Such  models  do  currently  exist 
and  will  be  discussed  in  a  later  section.  However,  most  hydrodynamic 
modelers  employ  a  spatial-averaging  technique  similar  to  the  turbulent 
time-averaging  technique  to  yield  either  one-  or  two-dimensional  models. 
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As  previously  noted,  the  present  state  cf  the  art  is  such  that  three- 
dimensional  models  are  not  practical  for  long-term  simulations  of  the 
hydrodynamics  of  a  reservoir. 

29.  The  basic  assumption  in  the  spatial  averaging  of  the  three- 
dimensional  equations  is  that  the  dependent  variables  can  ce  represented 
by  an  average  value  over  one  or  mere  of  the  spatial  coordinates  plus 
some  small  random  deviation;  e.g.,  the  velocity  would  be  written  as 

u .  =  u .  +  u !  (l3) 

111 

where 

x. +Ax. 1 2 
1  1 

u.  =  — —  /  u.  dx. 

1  Ax.  J  11 

1  x.-Ax./2 
1  1 

x.+Ax./2 
1  1 

/  u!  dx.  =  0 
Ax.  /  11 

1  x.-Ax./2 
1  1 

and 

u^  =  time-  and  space-averaged  velocity 
Ax^  =  spatial  step 

u!  =  deviation  between  time-averaged  velocity  and  time-  and 
space-averaged  velocity 

In  an  x  ,  y  ,  z  coordinate  system  (with  x  referring  to  the  longi¬ 
tudinal;  y  ,  the  lateral;  and  z  ,  the  vertical),  if  i  =  ?  ,  the  inte¬ 
gration  is  over  the  width  and  a  width-averaged  model  results.  However, 
if  i  =  3  ,  the  integration  is  taken  over  the  depth  and  a  depth-averaged 
model  will  result.  Many  depth-averaged  models  have  been  developed  since 
Leendertse's  (1967)  work;  whereas,  laterally  averaged  models  have  only 
been  developed  over  the  past  five  years  or  so.  If  the  integration  is 
performed  over  the  complete  cross  section,  a  one-dimensional  model  with 
variations  allowed  only  in  the  longitudinal  direction  results.  Such 
models  are  not  considered  in  this  study. 

30.  As  was  done  in  the  time-averaging  of  the  instantaneous 
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equations,  expressions  such  as  Equation  13  are  substituted  into  the 
turbulent  time-averaged  equations  to  yield  a  set  of  equations  with  the 
time-averaged  and  spatially  averaged  components  of  the  flow  and  thermo¬ 
dynamic  variables  as  des endent  variables  plus  the  additional  terms 


and 
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As  in  the  time-averaging  case,  these  terms  are  normally  approximated  by 
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where  e! .  and  A! .  are  referred  to  as  "eddy  dispersion  coefficients" 
ij  ij 

by  Holley  (1969)  to  distinguish  them  from  the  turbulent  eddy  diffusion 
coefficients  arising  from  the  time  averaging,  and  T  is  the  time- 
averaged  and  spatially  averaged  temperature. 

31.  The  resulting  spatially  averaged  equations  take  different 
forms,  depending  upon  whether  the  averaging  is  performed  over  the  depth 
or  the  width.  Since  depth-averaged  models  are  not  applicable  to  the 
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hydrodynamic  modeling  of  stratified  reservoirs,  only  the  laterally 
averaged  equations,  with  B  denoted  as  the  width,  are  presented  below: 


in  hydrodynamics,  is  to  write  the  equations  with  the  stream  function  and 
the  vorticity  as  the  dependent  variables.  With  the  vorticity  defined  as 
the  curl  of  the  velocity,  i.e.,  u  =  V  x  v  ,  it  can  easily  be  seen  that 
if  the  velocity  field  is  two-dimensional,  e.g.,  v  =  ui  +  vj  ,  where 
i  ,  j  and  ..  are  unit  vectors,  the  vorticity  contains  only  one  compo¬ 
nent,  namely  ck  .  Therefore,  instead  of  being  required  to  solve  two 
momentum  equations  for  the  velocities  u  and  v  ,  a  single  equation  for 
C  and  a  Poisson  equation  for  the  stream  function  are  solved.  In  other 
words,  the  number  of  equations  to  be  solved  has  in  essence  been  reduced 
by  one.  However,  one  must  make  additional  computations  to  obtain  the 
velocity  field  from  the  computed  stream  function.  Still,  when  the 
vorticity-stream  function  approach  is  applicable,  it  is  probably  faster. 
Multiple  outlets  at  a  dam  would,  however,  prohibit  its  use. 

Subgrid-scale  motion 

31*.  The  eddy  coefficients  discussed  above  enter  the  equations  due 
to  first  of  all  the  time  averaging  (diffusion  coefficients)  of  the  equa¬ 
tions  and  then  as  a  result  of  spatial  averaging  (dispersion  coefficients) 
to  remove  one  or  more  of  the  independent  spatial  coordinates  from  the 
equations.  A  similar  coefficient  arises  as  a  result  of  averaging  the 
governing  equations  over  the  numerical  grid  upon  which  a  numerical  solu¬ 
tion  is  sought.  The  numerical  model  cannot  resolve  small-scale  local 
circulation  patterns  or  eddies  unless  the  eddies  extend  over  an  area 
covering  several  grid  points.  Thus,  as  discussed  by  Deardorff  (1970), 
an  averaging  operator  is  applied  to  the  governing  equations,  with  aver¬ 
aging  typically  being  over  the  grid  volume  of  the  numerical  calculation 
to  filter  out  the  subgrid  scale  (SGS)  motions.  Explicit  calculations 
are  then  made  for  the  filtered  variables  after  assumptions  are  made 
about  the  SGS  Reynolds  stresses  that  arise  from  the  averaging  process. 

All  of  this,  of  course,  is  completely  analogous  to  the  manner  in  which 
turbulent  and  dispersive  Reynolds  stresses  arose  in  the  previously 
discussed  time  and  spatial  nv  -aging. 

35 •  The  total  stress  then  is  the  sum  of  the  molecular  viscous 
stress,  the  diffusive  Reynolds  stress,  the  dispersive  Reynolds  stress. 
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and  the  SGS  Reynolds  stress.  In  practice,  all  of  these  are  lumped  into 
one  stress  term  with  a  single  tensor  eddy  viscosity  coefficient.  Simi¬ 
larly,  the  total  diffusivity  in  the  transport  equation  for  temperature 
is  the  sum  of  four  components  that  are  lumped  together  with  a  single 
tensor  eddy  diffusivity  coefficient. 

3 6.  In  turbulent  flow,  these  coefficients  are  not  constant  as  in 
laminar  flow,  but  depend  on  the  flow  itself,  i.e.,  on  the  processes 
generating  the  turbulence.  The  determination  of  these  eddy  coefficients 
in  terns  of  the  mean  flew  variables  is  a  major  problem  in  hydrodynamic 
modeling. 

37.  Up  to  this  point,  the  eddy  viscosity  and  diffusivity  coeffi¬ 
cients  have  been  treated  as  second  order  tensors.  Some  researchers 
actually  allow  for  the  tensor  behavior  as  a  function  of  the  rate  of 
strain  tensor  of  the  mean  flow  field;  however,  the  more  common  approach 
is  to  neglect  all  off-diagonal  terms  and,  furthermore,  to  consider  the 
two  components  in  the  horizontal  plane  to  be  equal.  Some  modelers  allow 
for  a  variation  of  these  coefficients,  but  others  take  a  rather  simplis¬ 
tic  approach  and  treat  them  as  constants  over  the  computational  field. 

38.  As  noted  by  Lick  (1976),  the  vertical  eddy  coefficients  should 
vary  throughout  the  depth.  Causes  for  their  variations  are  related  to 
the  following: 

a.  Stability  of  the  water  column. 

b.  Action  of  the  wind  on  the  surface. 

c_.  Vertical  shear  in  currents  due  to  horizontal  pressure 
gradients . 

d_.  Presence  of  internal  waves. 

e_.  Effect  of  bottom  irregularities  and  friction  on  currents. 

One  often  finds  the  vertical  coefficients  related  to  the  stability  of 
the  water  column  as  a  function  of  the  Richardson  number  R.  : 


£  j£ 
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where  z  is  the  vertical  coordinate  and  u  is  the  mean  horizontal 
flow  velocity. 

39*  Horizontal  coefficients  are  generally  much  greater  than  verti¬ 
cal  coefficients.  Based  upon  various  surface  dye  experiments,  primarily 
in  the  oceans,  it  has  been  found  that  the  horizontal  coefficients  are 
proportional  to  the  scale  of  the  turbulence  raised  to  approximately  the 
U/3  power. 

1*0.  Lick  indicates  that  in  nonstratified  flow,  the  eddy  diffusiv- 
ity  is  approximately  equal  to  the  eddy  viscosity,  i.e.,  the  Reynolds 
analogy  holds.  Various  forms  that  have  been  employed  for  these  coeffi¬ 
cients  will  be  presented  later  in  discussions  on  individual  models. 

Boundary  Conditions 

1*1.  As  noted  by  Roache  (1972),  the  thing  that  makes  a  particular 
fluid  flow  problem  unique  are  the  boundary  conditions  that  are  prescribed. 
Conditions  at  the  surface  of  the  reservoir,  at  solid  boundaries,  and  at 
both  inflow  and  outflow  boundaries  must  be  specified  in  order  to  obtain 
a  solution  of  the  governing  equations  previously  presented. 

Surface  conditions 

1*2.  In  modeling  the  hydrodynamics  of  a  water  body,  one  of  two 
approaches  is  taken  in  the  treatment  of  the  water  surface.  The  surface 
is  either  treated  as  a  free  surface  or  as  a  rigid  lid.  In  either  case, 
the  heat  flux  at  the  surface  must  be  specified  as  a  boundary  condition 
on  the  temperature. 

1*3.  If  the  surface  is  treated  as  free,  the  assumption  is  made 
that  a  water  particle  on  the  surface  remains  there,  i.e.,  the  surface 
is  a  streamline.  This  then  leads  to  the  following  kinematic  condition: 

f+UU+Vf-W-°  (l8) 

for  the  computation  of  the  water  surface  elevation,  C  .  In  addition, 
the  internal  stresses  in  the  fluid  must  equal  those  applied  externally 
to  the  surface.  Considering  a  vertical-longitudinal  two-dimensional  (2-D) 


reservoir,  for  laminar  flow  the  normal  internal  stress  at  the  surface 


Figure  3.  Orientation  of  unit  normal  to  the  surface 


components,  respectively,  of  the  outward  unit  normal  vector  to  the  sur¬ 
face.  The  above  expressions  have  been  derived  from  the  stress  force 
t  given  by 


where,  as  noted,  n  is  the  unit  normal  to  the  surface  and  T  is  the 
laminar  stress  tensor  for  incompressible  flow,  given  by 


Now  the  externally  applied  stress  will  be  a  normal  stress  as  a  result  of 
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atmospheric  pressure  and  a  tangential  stress  imparted  by  the  wind; 
therefore,  the  boundary  conditions  become 


and 
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where 

T, =  wind  shear  stress 

iv  i  r< ;  i 

P&  =  atmospheric  pressure 

Thus,  in  a  strict  application  of  the  stress  boundary  condition,  the 

orientation  of  the  surface,  i.e.,  n  and  n  ,  would  have  to  be  known. 

x  z 

Since  for  a  large  water  body  the  surface  is  at  least  locally  flat,  the 
assumption  of  a  flat  surface  is  normally  made  so  that 
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The  stress  boundary  conditions  then  reduce  to 


P  =  P  - 
a 


UU.  In  addition,  if  the  hydrostatic  pressure  assumption  is  made, 
i.e.,  vertical  accelerations  are  neglected,  the  above  conditions  take 
the  form  below  that  is  commonly  found  in  the  literature,  where  the 
molecular  viscosity  y  has  been  replaced  by  its  turbulent  eddy 

counterpart ,  . 
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U5.  When  the  surface  is  treated  as  a  rigid  lid,  it  becomes  in 


essence  a  solid  boundary,  and  the  normal  component  of  the  velocity  must 
be  zero.  In  addition,  the  pressure  can  no  longer  be  prescrined  at  the 
surface,  but  rather  must  be  computed.  The  pressure  boundary  condition 
now  takes  the  form  of  a  derivative  boundary  condition,  i.c-.,  a  Neummat. 
condition  as  opposed  to  the  Dirichlet  condition  for  the  free  surface 
case.  The  primary  reason  for  assuming  that  the  surface  is  a  rigid  lid 
is  connected  with  stability  problems  encountered  in  the  numerical  solu¬ 
tion  of  the  governing  equations  and  will  be  discussed  later.  Tt  should 
be  obvious  that  with  the  rigid  lid  approximation,  the  effect  on  the 
internal  flow  of  the  piling  up  of  water  cannot  be  accounted  for. 

Modelers  such  as  Liggett  (1970)  and  Lick  (1976)  have  employed  the  concept 
in  the  development  of  lake  circulation  models.  However,  others  such  as 
Eraslan*  and  Edinger**  feel  that  the  rigid  lid  approximation  associated 
with  a  uniform  water  surface  assumption  is  not  realistic  in  the  develop¬ 
ment  of  mathematical  models  for  environmental  flow  conditions  and  that 
the  water  surface  elevation  must  be  considered  as  an  integral  part  of 
the  general  solution  of  the  hydrodynamic  problem. 

Solid  boundaries 

U6.  For  viscous  fluids,  the  fluid  velocity  is  actually  always 
zero  at  a  solid  boundary;  i.e.,  both  the  tangential  as  well  as  the  ncimal 
components  are  zero.  This  boundary  condition  is  referred  to  as  a  "nc- 
slip  condition."  Although  in  theory  such  a  condition  must  always  hold 
at  a  solid  boundary,  often  in  hydrodynamic  modeling  a  slip  condition  is 
employed.  This  condition  is  implemented  by  setting  the  component  of  the 
velocity  normal  to  the  wall  equal  to  zero  but  not  the  tangential;  i.e., 
the  flow  slides  freely  along  the  solid  wall.  Theoretically,  this  implies 
that  the  flow  is  inviscid.  Therefore,  proper  boundary  conditions  for 
a  slip  wall  are  that  the  normal  velocity  is  zero  as  well  as  that  the 
vorticity  is  zero  at  the  solid  wall,  since  vorticity  is  created  in 


*  Personal  communication.  May  1979,  Arsev  Eraslan,  Chief  Scientist, 
Hennington,  Durham,  and  Richardson,  Knoxville,  Tenn. 

**  Personal  communication.  May  1979,  J.  E.  Edinger,  J.  E.  Edinger 
Associates,  Inc.,  Wayne,  Penn. 
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viscous  regions.  The  tangential  component  cf  velocity  is  then  deter¬ 
mined  from  ttie  zero-vort ic i ty  condition.  However,  it  does  net  appear 
that  this  is  the  approach  usually  taken  ir.  hydrodynamic  modeling.  In 
most  crises,  as  will  he  discussed  later,  a  staggered  grid  is  used  in 
obtaining  a  numerical  solution  such  that  the  tangential  component  of  the 
velocity  is  not  defined  at  the  wall.  Instead,  its  value  must  he  speci¬ 
fied  outside  the  wall.  The  usual  procedure  taken  fcv  most  hydrodynamic 
modelers  for  slip  walls  is  to  set,  this  value  equal  to  its  value  inside 
the  field. 

1*7.  The  major  reason  for  using  clip  boundary  conditions  is 
apparently  related  to  the  fact  that  a  relatively  large  grid  spacing  is 
normally  required  in  hydrodynamic  modeling  for  economic  reasons. 
such  a  grid  spacing  near  a  solid  boundary,  if  nc-sl ’p  conditions  are 
applied,  t  he  boundary  layer  effect  extends  farther  into  the  field  t  ‘.an 
it  does  in  reality. 

I48.  In  addition  to  conditions  being  imposed  on  the  flew  field  at 
solid  boundaries,  information  about  the  heat  transfer  must  also  le 
specified.  Either  wall  temperatures  or  the  heat  flux  may  be  preserved. 
In  all  reservoir-  or  lake-type  modeling  that  has  been  investigated,  the 
solid  boundaries  are  assumed  to  be  adiabatic,  and  thus  the  heat  flux 
through  the  boundary  is  set  tc  zero. 

Open  boundaries 

1*9.  Open  boundaries  are  exactly  what  the  name  implies,  i.e., 
boundaries  that  are  open  such  that  fluid  may  either  enter  or  leave  the 
field  within  which  a  solution  is  sought.  Tuch  boundaries  are  known  as 
either  "inflow"  or  "outflow"  boundaries. 

90.  At  forced  open  boundaries  (inflows  are  always  forced),  cither 
flow,  i.e.,  velocities,  or  water  surface  elevations  (assuming  a  free  sur¬ 
face),  must  be  prescribed  as  a  function  of  time  along  with  the  tempera¬ 
ture.  Theoretically,  rather  than  expressing  either  the  flow  or  surface 
elevations,  one  could  specify  a  relationship  between  the  two.  fuch  a 
boundary  condition,  known  as  a  rating  curve,  is  often  prescribed  at  the 
downstream  end  of  one-dimensional  unsteady  flow  models  of  river  flows. 

51.  At  outflow  boundaries  that  are  free  rather  than  forced,  e.g.. 
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free  flow  from  an  opening  at  a  dam,  one  cannot  prescribe  the  flow 
directly  since  it  obviously  is  dependent  upon  conditions  within  the 
computational  field.  In  the  temperature  computations,  an  outflow  bound 
ary  is  always  considered  free.  A  boundary  condition  that  allows  phenom 
ena  generated  in  the  domain  of  interest  to  pass  through  the  boundary 
without  undergoing  significant  distortion  and  without  influencing  the 
interior  solution  is  needed.  Since  a  physical  law  to  prescribe  such  a 
boundary  condition  does  not  exist,  some  kind  of  extrapolation  from  the 
interior  must  be  used.  The  most  common  methods  used  are  either  a 
Sommerfeld  radiation  condition  or  perhaps  one-sided  differences  when 
employing  finite  differences  to  obtain  numerical  solutions. 

r>?.  The  dispersion  characteristics  in  one  dimension  of  the  waves 
needed  to  prescribe  the  Commerfeld  radiation  condition  are  known  as 
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where  $  is  any  variable  and  C  is  the  :  velocity  f  the  waves. 

The  dispersion  characteristics  are  r.-'t  general  iy  known,  since  '  is 
not  generally  a  constant.  A  simpli  f  imt  ion  'ft  er  as>-d  *  nv  j  ‘his 
i  rebletn  is  to  assume  that  the  ohara  “••!•  i  si  i  -s  •'  ‘is  waves  will  huv-  a 
slope  e  )ual  to  Ax /At  (Ax  and  t  ro  I  a  i  -in  i  • .  m:  ••a. 

siet.  ,  r--S;  ect  i  vely  )  .  Then  ti.e  •  xtra;  '.•!’••}  var  •  :  a*  t  • 

•  irr.e  ‘  at  the  boundary  :  oin*  is  c . 


A  more  aeeurate  use  f  ‘tie  .'omtier  *'•  ■ .  !  ri  i  ;•<*  ti‘  i  :  .a*  i  \ 

is.  presented  by  ri.ans.k:  ,  l't.'t  .  las',  at  .  f  ,-i  >:  it.,-  a  ■  :.s*  a:.‘  v.  i  . 
for  ’tie  phase  velocity  ,  rianski  nurc-r  i  i  ...  a.  1  a*  es  a  -  r  ;  ura-  i 
velocity  from  neighboring 
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PART  III:  NUMERICAL  DISCUSSIONS 


53-  The  governing  equations  of  fluid  motion  are  nonlinear  partial 
differential  equations,  which  in  a  strict  mathematical  sense,  are  classi¬ 
fied  as  being  of  the  parabolic  type.  However,  outside  the  boundary 
layer  the  equations  exhibit  a  strong  hyperbolic  or  wave  character  due  to 
the  convective  terms  and,  thus,  are  often  considered  as  being  of  the 
hyperbolic  type.  In  any  case,  because  of  their  nonlinearity,  analytical 
olutions  do  not  generally  exist  and  one  must  resort  to  numerical  methods 
to  obtain  an  approximation  of  the  continuous  solution  of  the  differential 
equations.  Such  methods  include  the  use  of  finite  differences  and  finite 
elements . 


Finite  Element  Method  (FEM) 


5^.  In  the  finite  element  approach,  the  field  is  divided  into 
finite  elements,  and  the  solution  is  approximated  by  a  chosen  function  on 
each  element.  This  function  contains  free  parameters,  which  are  evalu¬ 
ated  by  requiring  the  function  and  perhaps  certain  of  it s  derivatives  to 
equal  the  solution  and  its  derivatives  at  certain  points  on  the  element, 
if  the  partial  differential  equations  can  be  expressed  in  terms  if  inte¬ 
gral  variational  principles,  tne  variational  integrals  over  eucn  element 
are  evaluated  analytically  from  tne  chosen  approximation  fund  Ions  :m 
each  element.  The  integrals  ver  each  individual  element  are  t  .hen  summed 
over  all  tne  e  iements  t  p  roduce  the  variational  integral  over  tne  i  re 
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functions  and  the  partial  differential  equations  are  evaluated  on  each 
element.  This  produces  a  set  of  simultaneous  algebraic  equation  to  be 
solved  for  the  values  of  the  solution  and  perhaps  some  of  its  derivatives 
at  certain  points  on  the  elements. 

55-  The  finite  element  approach  is  best  suited  to  partial  differ¬ 
ential  systems  that  can  be  expressed  in  terms  of  a  variational  principle. 
In  this  case,  the  boundary  conditions  can  be  incorporated  naturally  by 
means  of  Lagrange  multipliers.  For  more  general  systems,  particularly 
nonlinear  systems  that  are  not  expressible  in  terms  of  variational  prin¬ 
ciples,  the  finite  element  approach  must  use  the  method  of  weighted 
residuals  (dalerkin)  whereby  a  functional  form  of  the  solution  in  each 
element  is  assumed  and  integral  moments  of  the  partial  differential  equa¬ 
tions  are  satisfied  over  the  field  as  noted  above.  With  this  procedure, 
the  partial  differential  equations  themselves  are  not  actually  satisfied. 
Boundary  conditions  are  incorporated  in  the  assumed  functional  form  of 
the  solution  in  the  elements  adjacent  to  the  boundaries. 

56.  The  finite  element  method  has  enjoyed  its  greatest  success  in 

the  field  of  solid  mechanics  where  for  the  most  part  variational  rather 

than  difference  met ho  is  are  used.  As  Fix  (1175)  notes : 

"The  reason  for  this  is  partly  physical.  The  equa¬ 
tions  of  elasticity  can  be  put  into  a  variational 
form  and  engineers  nav  ■  found  this  to  be  the  most 
physically  natural  setting  to  formulate  apt  roximat ions . 

In  addition,  the  variational  approximations — finite 
elements — have  other  properties  that  are  of  great, 
value  in  practice.  Complicated  boundaries  can  easily 
:  ••  treated  in  this  setting;  singui unifies  in  the  solu¬ 
tion  can  tie  modeled  in  tne  approximation;  and,  in 
leading  wifi,  higher  order  met nous  vltn  increased  re- 
ivinc  ;  .vr,  the  practical  problems  are  much  less 
t  r . 'Utleso::*-  than  wi  si  difference  schemes." 

1' .  :  er.nuj  .  ’:,e  major  practical  disadvantage  of  the  finite  element 

met  :.o  j  as  'ijp.ei  to  hydrodynamic  p  roblems  1.,  that  a  large  amount  of  eom- 
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matrices  that  result  from  using  differences.  A  common  rebuttal  by  finite 
element  modelers  to  this  criticism  is  that  since  implicitness  is  central 
to  the  finite  element  approximation,  implicit  time  differencing  should 
be  employed  to  yield  an  unconditionally  stable  system  and  thus,  the  use 
of  larger  time  steps  compensates  for  the  large  computing  times  required 
to  invert  the  dense  coefficient  matrices  each  time  step.  However,  it 
should  be  remembered  that  in  addition  to  stability  considerations,  one 
must  first  of  all  be  concerned  with  the  accuracy  of  the  solution.  When 
using  finite  differences,  it  can  be  shown  that  as  the  computational  time 
step  becomes  increasing  larger  than  that  allowed  by  the  Courant  condi¬ 
tion  for  gravity  waves,  e.g.,  in  one  dimension 

At  < 

C  ^h 


where 

At^  =  time  step  restricted  by  Courant  condition 
h  =  water  depth 

a  corresponding  increase  in  the  number  of  spatial  points  per  wavelength 
must  occur  to  retain  the  same  level  of  accuracy  in  the  amplitude  and 
phase  of  the  computed  wave.  Leendertse  (1967)  indicates  that  from  a 
practical  standpoint,  generally  the  time  step  should  not  be  greater  than 
3-5  times  At^  in  the  difference  scheme  he  employs  to  compute  vertically 
averaged  flows.  Abbott  (1979)  and  his  colleagues  at  the  Danish  Hydraulic 
Institute  (DHI)  have  arrived  at  a  similar  conclusion  for  the  difference 
scheme  that  they  employ.  Edinger  and  Buchak  (1979)  have  indicated  that 
for  the  laterally  averaged  case,  the  time  step  should  not  exceed  3-5 
times  the  At  computed  from  the  internal  wave  condition  given  by 


At  < 


V 


Ax _ 

Ap  gh 
P 


where 

Ap  =  change  in  water  density 
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Such  an  analysis  of  the  error  between  the  computed  and  analytic  wave  in 
the  FEM  has  not  been  found  in  the  literature. 

58.  In  the  de’ elopment  of  a  new  model,  an  often  stated  disadvan¬ 
tage  of  the  finite  element  method  is  that  the  FORTRAN  coding  is  much 
more  cumbersome  than  with  finite  differences.  Of  course,  if  the  model 
already  exists,  such  a  disadvantage  is  of  no  concern  to  the  user  unless 
extensive  modifications  are  required,  in  which  case  the  cumberscmeness 
of  the  coding  might  well  become  a  major  consideration  in  the  model 
selection. 


Finite  Difference  Method  (FDM) 


59-  The  vast  majority  of  the  numerical  hydrodynamic  models, 
whether  they  be  one-,  two-,  or  three-dimensional,  employ  the  use  of 
finite  differences  to  obtain  solutions  of  the  governing  equations  of 
fluid  motion.  In  the  finite  difference  method,  the  domain  of  the  indet 
dent  variables  is  replaced  by  a  finite  set  of  points  referred  to  as  net 
or  mesh  points.  One  then  seeks  to  determine  approximate  values  for  the 
desired  solutions  at  these  points.  The  values  at  the  mesh  points  are 
required  to  satisfy  difference  equations  that  are  usually  obtained  by 
replacing  partial  derivatives  by  partial  difference  quotients.  The  re¬ 
sulting  set  of  simultaneous  algebraic  equations  is  then  solved  for  tne 
values  of  the  solution  at  the  mesh  points.  If  the  boundaries  do  not 
coincide  with  mesh  points,  then  the  finite  difference  approach  applied 
to  the  equations  in  a  Cartesian  coordinate  system  requires  interpolatio 
between  mesh  points  to  represent  boundary  conditions. 

bO.  However,  through  coordinate  transformations ,  irregular  bound 
aries  can  be  accurately  handled  while  still  making  use  of  the  sir.pl  icit 
of  finite  differences  to  obtain  solutions.  The-  most  general  of  rue:, 
transformations,  which  will  he  discussed  in  more  detail  Inter  in  tin  ’-e 
port,  is  a  method  developed  by  Thompson  et  al.  which  nenorates 

curvilinear  coordinates  as  the  solution  of  two  elliptic  partial  differs 
tial  equations  with  Dirichlet  boundary  conditions,  one  coordinate  reins 
specified  as  constant  on  the  boundaries,  and  a  distribution  of  the  el.no 
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specified  along  the  boundaries.  No  restrictions  are  placed  on  the  irreg¬ 
ularity  of  the  boundaries,  and  fields  containing  multiple  bodies  or 
branches  can  be  handled  as  easily  as  simple  geometries.  Regardless  of 
the  shape  and  number  of  bodies  and  regardless  of  the  spacing  of  coordi¬ 
nate  .  .nes,  all  numerical  computations,  both  to  generate  the  coordinate 
system  and  +  j  subsequently  solve  the  fluid  flow  equations,  are  done  on 
a  rectangular  grid  with  square  mesh. 

61.  Since  the  boundary-fitted  coordinate  system  has  a  coordinate 
line  coincident  with  all  boundaries,  all  boundary  conditions  may  be 
expressed  at  grid  points,  and  normal  derivatives  may  be  represented 
using  only  finite  differences  between  grid  points  on  coordinate  lines. 

No  interpolation  is  needed,  even  though  the  coordinate  system  is  not 
orthogonal  at  the  boundary. 

62.  Linear  transformations  that  allow  for  the  physical  dimensions 
to  be  mapped  between  the  values  of  0  and  1  have  been  employed.  For 
example,  as  will  be  discussed  later  in  PART  IV,  Lick  (1976)  maps  the 
vertical  dimension  in  such  a  manner  to  represent  bottom  topographies 
more  accurately. 

Discrete  element  concept 

63.  Eraslan  employs  a  numerical  technique  that  he  labels  "the 
discrete  element  method."  However  from  a  conceptual  standpoint,  the 
primary  difference  between  the  finite  difference  method  as  it  is  normally 
applied  and  the  discrete  element  method  appears  to  be  that  the  mathemati¬ 
cal  development  of  the  discrete  element  method  is  based  on  employing  the 
control  volume  integral  forms  of  the  physical  conservation  principles; 
whereas,  the  usual  application  of  the  finite  difference  method  begins 
with  the  continuum  limit  differential  equations  presented  in  PART  II. 

6L.  Eraslan  indicates  that  the  application  of  the  discrete  element 
method  to  the  solution  of  environmental  fluid  mechanics  problems  is  based 
on  the  following  procedure: 

a.  Divide  the  flow  region  into  arbitrarily  sized  discrete 

elements,  preferably  with  geometrically  simple  (rectangu¬ 
lar)  enclosure  surfaces  except  at  the  boundaries,  such 
that  the  finite  number  of  discrete  elements  completely 
spans  the  region. 
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b. 


Integrate  the  volume  and  surface  area  integrals  of  the 
physical  conservation  equations  without  assuming  uniform 
values  for  the  flow  properties  over  the  surface  areas. 
This  produces  a  governing  serai  discretized  system  of  ordi¬ 
nary  differential  equations  in  time. 

e_.  Apply  proper  interpolation  techniques  for  determining 
transportive  values  of  trie  flow  properties  between  dis¬ 
crete  elements. 

6p.  As  an  example  of  the  discrete  element  concept,  consider  the 
one-dimensional  problem  illustrated  in  Figure  lr.  ,'ieglecting  frictional 


Figure  4.  One-dimensional  discrete  element 


effects,  the  integral  forms  of  the  conservation  of  mass  and  momentum  can 
be  written  as 


Continuitj 


v  •  n  dA  =  0 


Momentum:  — 

O  t, 


Iff  *  «.  *  ffA  ” 

J  J  J  cv  J  cv 


n  dA  = 
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where 


dV  =  differential  volume 
o 

cv  =  volume  of  discrete  element 

A  =  Area  of  discrete  element 
cv 

> 

v  =  velocity  vector 

clA  =  differential  area 
> 

f  =  force  vector 


Now  define 


u  dA  =  volumetric  flow  rate 
A 

cv 


Therefore,  from  Figure  U 


IL 


vndA=  -Gi_1/2  +  Gi+l/2 


cv 


also. 


3 


Iff 


3H. 


n ,  111  dV  =  f-  (Ax. A.)  =  Ax . B .  ~ 

3t  III  o  3t  1  1  l  l  at 


cv 


where  B.  is  the  surface  width  and  H.  is  the  surface  elevation. 

l  i 

Therefore,  the  discrete  element  equation  for  the  conservation  of  fluid 
mass  becomes 


3t  =  AxiBi  (Gi-l/2  °i+l/2^  ^ 

Considering  the  surface  integral  for  the  momentum  flux  over  the  cross 
section  A^  at  the  center  of  the  element  yields 

vv  •  n  dA  =  +  Reynolds  stress  terms 
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Neglecting  the  Reynolds  stress  terms,  which  arise  due  to  writing  the 
velocity  as  the  sum  of  a  time-averaged  and  spatially  averaged  component 
plus  its  time  and  spatial  deviations,  the  integral  form  of  the  momentum 
equation  applied  to  one-half  of  the  discrete  element  becomes 


66.  Equations  22  and  23  both  take  the  appearance  of  finite  differ¬ 
ence  equations  in  which  the  time  derivative  has  not  been  replaced  by 
differences.  It  appears  that  from  a  practical  consideration,  the  primary 
difference  between  the  discrete  element  method  and  the  application  of 
finite  differences  to  the  differential  equation  centers  around  what  might 
be  called  "the  conservation  of  geometrical  properties,"  as  reflected 
through  the  definition  of  the  divergence  of  a  variable.  In  the  equations 
of  motion,  flux  terms  such  as 


f  f  $  •  n  dA 

J  J  ^cv 


appear.  Considering  the  element  below. 


and  working  with  only  the  x  direction,  the  flux  integral  above  can  be 
evaluated  as 


ILJ  • : 


dA  =  (4>A)+  -  (<*>A)_ 


Now,  if  one  employs  Gauss's  Divergence  theorem,  the  flux  integral  can 
be  written  as 


/ / »ov5 ' " “ '  IIL  (7  •  $)  dVQ 


s  (?  *  4')idxiAi 


where  7  •  $  is  the  divergence  of  $  .  The  two  previous  expressions 
can  be  equated,  and  one  can  derive  an  expression  for  the  divergence  over 
the  i^h  element  as 


('  •  *)i  '  ITT  ><♦*>♦  -  <«>J 


('  •  Oi  -  ^ 


In  the  usual  derivation  of  the  differential  form  of  the  equations,  the 
divergence  is  written  as 

y  .  *  —  (pc 


Equation  24  might  be  referred  to  as  the  "geometrically  conservative" 
form  of  the  divergence;  whereas,  the  normal  definition  as  given  by  Equa¬ 
tion  25  would  be  referred  to  as  the  "geometrically  nonconservative"  form. 

6j.  Note  that  physically  the  difference  between  the  two  forms  is 
that  in  the  conservative  form  (Equation  24),  the  area  through  which  the 
flux  of  flows  is  that  of  the  bounding  faces  through  which  the  flux 
actually  occurs.  In  Equation  25,  however,  the  influence  of  the  respec¬ 
tive  areas  on  the  flux  through  the  boundaries  is  not  allowed. 

68.  It  appears  that  if  the  conservative  form  of  the  divergence  is 


used  in  expanding  the  vector  form  of  the  governing  dif ferort ial  equa¬ 
tions,  the  finite  difference  method  applied  to  those  equations  becomes 
identical  to  the  discrete  element  method  if  the  same  interpolation  scheme 
is  used  to  provide  values  of  the  dependent  variables  at  points  where 
they  are  not  defined  in  the  grid. 

Finite  difference  spatial  grids 

69.  The  spatial  grid  in  Cartesian  coordinates  most  commonly  used 
by  numerical  hydrodynamic  modelers  appears  to  be  one  in  which  the  water 
surface  elevation,  temperature,  and  density  are  defined  at  the  center  of 
a  computational  cell;  whereas,  the  velocity  components  are  defined  on 
the  faces  of  the  cell.  Such  a  grid  is  illustrated  below  for  a  two- 
dimensional  problem. 


h,  p,  T 


With  such  a  grid,  the  normal  component  of  the  velocity  at  solid  bounda¬ 
ries  can  easily  be  set  to  zero  if  the  boundary  is  assumed  to  lie  along 
cell  faces,  which  is  the  usual  assumption. 

TO.  With  such  a  grid  one  obviously  will  need  values  of  variables 
at  points  where  they  are  not  defined  in  order  to  nmerically  solve  the 
governing  equations.  One  solution  is  to  utilize  more  than  one  grid, 
with  the  variables  defined  such  that  a  solution  on  one  grid  is  used  to 
step  the  solution  forward  on  another  grid,  fis  discussed  by  Simons 
(1973),  such  a  procedure  can  result  in  serai -independent  solutions  on  the 
different  grids.  The  numerical  error  associated  with  the  use  of  more 
than  one  grid  is  known  as  a  "grid  dispersion  error."  The  approach 


normally  taken  to  provide  the  variables  at  net  points  where  they  are  not 
computed  is  to  perform  an  interpolation  within  the  grid.  Simons  shows 
that  conservation  requirements  are  satisfied  if  the  unknown  values  are 
approximated  by  simple  linear  interpolation. 

71.  A  grid  often  used  in  aerodynamic  flow  modeling  has  all  vari¬ 
ables  defined  at  the  same  point,  i.e.,  at  the  cell  center.  Such  a  grid 
(shown  below)  has  been  employed  by  Waldrop  and  Tatom  (1976)  in  their 
three-dimensional  hydrodynamic  modeling  work. 


72.  Still  another  grid  is  currently  being  employed  by  Thompson* 
in  the  development  of  a  model  for  use  in  selective  withdrawal  studies. 


J  +  1/2  — — , 

u,  V 

J  -  1/2  — i 

1 

1-1/2  I  I  +  1/2 


*  Personal  communication,  April  1979,  J.  F.  Thompson,  Mississippi  State 
University,  Mississippi  State,  Miss. 


In  this  grid,  all  velocity  components  are  defined  at  the  same  point, 
i.e.,  the  cell  corners;  whereas,  all  thermodynamic  variables  are  defined 
at  the  cell  center.  Thompson  indicates  such  a  grid  allows  for  a  more 
natural  application  of  velocity  and  pressure  boundary  conditions  in  a 
curvilinear  coordinate  system. 

73*  Most  numerical  finite  difference  hydrodynamic  models  employ  a 
constant  grid  size  in  each  direction.  However,  models  have  been  devel¬ 
oped  that  allow  for  the  size  of  the  computational  cell  to  vary  over  the 
region  within  which  flow  computations  are  being  made  in  order  to  increase 
the  resolution  in  certain  areas.  Examples  are  the  3-D  models  of  Tatom 
and  Waldrop  and  Thompson's  2-D  model  that  utilizes  boundary-fitted  coor¬ 
dinates.  As  discussed  by  Roache  (1972),  there  are  two  approaches  to  the 
implementation  of  a  variable  computational  mesh.  One  can  merely  solve 
the  given  equations  on  a  grid  that  has  physically  been  constructed  such 
that  the  computational  nodes  are  not  evenly  spaced,  or  one  can  transform 
the  equations  and  solve  them  in  a  transformed  rectangular  plane  with 
equal  grid  spacing,  although  the  grid  spacing  is  not  equal  over  the  phys¬ 
ical  region.  Even  though  the  two  approaches  might  appear  to  be  similar, 
Roache  indicates  they  are  fundamentally  different.  When  the  untrans¬ 
formed  equations  are  differenced  in  the  variable  mesh,  the  result  is  a 
deterioration  of  formal  accuracy,  but  the  transformed  equations  may  be 
differenced  in  a  regular  mesh  with  no  deterioriation  in  the  formal  order 
of  truncation  error  relative  to  the  transformed  plane.  Roache,  there¬ 
fore,  states  that  the  coordinate  transformation  approach,  which  can  be 
used  for  the  purpose  of  aligning  coordinates  along  physical  boundaries 
as  well  as  increasing  resolution  in  certain  areas,  is  to  be  preferred. 

As  previously  noted,  Thompson's  boundary-fitted  coordinate  technique 
provides  the  most  general  such  transformation  that  can  be  attained. 


Time  differencing 


7*+.  As  previously  noted,  time  integration  is  performed  by  finite 


differences  even  in  the  finite  element  method,  ouch  time  differencing 


can  basically  be  classified  as  either  explicit  or  implicit.  For  either 


type,  one  can  construct  first,  second,  or  even  higher  order  schemes; 
although  Kreiss  (1975)  indicates  that  second  order  schemes  are  to  be 


preferred  if  one  combines  accuracy  with  considerations  of  economy  and 
simplicity.  As  an  additional  classification,  one  often  finds  time- 
differencing  schemes  referred  to  as  one-  or  two-step  schemes. 

75-  In  order  that  these  concepts  may  be  better  understood,  con¬ 
sider  the  following  basic  equation: 


+ 

at 


(26) 


76.  If  u  is  zero,  this  equation  is  the  parabolic  time-dependent 
diffusion  equation  in  which  the  dependent  variable  $  can  change  only 
through  the  second  order  derivative  dissipative  term.  If  the  diffusion 
coefficient  a  is  zero,  the  equation  is  a  hyperbolic  wave-type  equation 
in  which  <J>  can  vary  only  through  advection  by  the  velocity  u  . 

77-  Assuming  that  <j>  is  continuous  and  possesses  continuous 
derivatives,  a  Taylor  series  expansion  in  time  yields 

2  2 

<t>(t  +  At)  =  <j>(t)  +  At  |f  +  ^5—  +  0(At3)  (27) 


thus,  one  can  solve  for  3<t>/3t  as 


ii  =  lit  1  At)  -i(tj  +  At  afi  +  0{At2} 

3t  At  2  « ,  2 

ox 


or 


If  =  +  0(At)  (28) 

dt  At 

This  is  called  a  "forward  difference"  representation  of  the  time  deriva¬ 
tive  and  as  indicated  by  the  notation  O(At)  is  only  first  order. 
Likewise,  one  can  write  the  Taylor  series  as 

2  2 

*{t  -  At)  =  4>(t)  "  |f  +  ff-  ^-f  +  0(At3)  (29) 

9t 
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so  that 


df  _  f  t  ;  -  ji  (  t  -  At, 
it  ”  At 


1 ;  At  ) 


which  id  Known  as  a  "backward  difference."  Ac-  with  the  forward  differ¬ 
ence,  such  an  expression  is  only  a  first  order  so hen.e .  If  one  subtracts 
Equation  29  front  Equation  27,  the  following  resuitr 

^  +  At)  -  tit  -  At)  [  ,  M  /  u ) 

3t  2At  U-  ' 

This  expression  is  referred  to  as  a  "centered  difference"  representation 
and  is  a  more  accurate  integration  scheme  as  At-n)  ,  since  .it  is  of 
second  order  in  time. 

78.  Applying  a  forward  differencing  of  the  time  derivative  in 
Equation  26  yields 


<fg  +  At 


1£  + 
3x 


a 


(32) 


which  is  Known  as  an  "explicit  time-integration  scheme,"  since  values  at 
the  n  +  1  time  level  can  be  computed  directly  from  known  values  at 
the  previous  time  level  n  .  In  addition,  such  a  scheme  is  labeled  as 
a  one-step  scheme,  since  only  one  sequence  of  computations  is  required. 

79-  Likewise,  applying  a  backward  differencing  to  the  time 
derivative  yields 


+  At 


ii  + 

3x 


3  2<t> 
3x2 
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(33) 


which  is  labeled  as  an  "implicit  time-integration  scheme,"  since  values 
at  the  ith  spatial  point  on  the  n  +  1  time  level  are  dependent  upon 
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(37) 


i£  + 

3t 


c  11  = 
3x 


0 


As  previously  illustrated,  a  centered  time  differencing  yields 


=  4>n~>'1  -  4>n  1  At2 
3t  -  At  3  at3  + 


0(At3) 


whe  five  time  levels  would  be  required  for  a  centered  difference  repre¬ 
sentation  of  the  third  order  time  derivative.  However,  through  a  sequence 
of  time  and  spatial  differentiation  of  Equation  37,  one  can  show  that 


lit  =  c3  ill 

3  °  3 

3tJ  3xJ 


thus,  a  third  order  time-differencing  scheme  of  the  basic  equation  could 
be  written  as 


q  ?  t 

d±  CJAt  2_i 
.  ^  3 

The  only  numerical  hydrodynamic  model  found  in  the  literature  that  is 
of  higher  order  than  two  in  time  is  a  vertically  averaged  model  devel¬ 
oped  at  the  Danish  Hydraulic  Institute,  which,  according  to  Abbott  (1979), 
is  "close"  to  third  order. 

Space  differencing 

83-  As  in  the  discussion  on  time  differencing,  either  first, 
second,  or  higher  order  differencing  of  spatial  derivatives  can  be  uti¬ 
lized  to  create  different  order  finite  difference  schemes.  Once  again, 
Taylor  series  expansions  in  space  yield  the  following  expressions  for 
forward,  backward,  and  centered  differences  of  a  spatial  derivative: 


1*9 


backward : 


,  3<t»  _  'I*i+1  * i-i  .  .,/  .  c 
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Ihus,  one  can  see  that  if  a  forward  difference  is  useu  in  both  space 
and  tine,  a  scheme  that  is  completely  first  order,  i.-.,  0;At,ux)  , 
results.  Likewise,  the  use  of  centered  differences  in  botu  space  and 
tine  results  in  an  0(At^,Ax^)  scheme. 

84 .  As  with  centered  differencing  of  tine  derivatives,  the  use  oi 
centered  difference!:  to  replace  spatial  derivatives  can  result  in  a  com¬ 
putational  node.  This  is  illustrated  by  considering  the  solution  of 


f  »  0 

9x 


where  the  analytic  solution  is  f  =  constant  .  Thus,  similar  to  its 
time  counterpart,  numerically  the  following  solution  can  develop: 


2  /'  ''  /\ 
f  /  w  w  y  \ 


POSSIBLE  NUMERICAL  SOLUTION 


■ANALYTIC  SOLUTION 


i  -  1  i  i  +  1 


85.  A  major  problem  associated  with  space  differencing  is  the 
treatment  of  the  nonlinear  advective  terms.  The  nonlinear  terms  tend  to 
generate  higher  harmonics,  which  can  result  in  what  Phillips  (19^9)  called 
a  "nonlinear  computational  instability."  As  noted  by  Roache  (1972), 
this  problem  is  not  unique  to  nonlinear  systems,  but  can  occur  whenever 
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nonconstant  coefficients  appear  in  the  differential  equation.  An  effec¬ 
tive  method  to  suppress  higher  harmonics  is  to  introduce  eddy  diffusion 
or  smoothing.  One  scheme  that  tends  to  introduce  artificial  damping  of 
the  higher  harmonics  without  appreciably  affecting  the  long  waves,  i.e., 
the  solution  of  interest,  is  the  two-step  Lax-V/endrof f  scheme  which  com¬ 
bines  Lax's  forward-in-time  and  eentered-in-space  scheme  as  the  first 
step  with  a  eentered-in-time  and  c enter ed- in-space  scheme  for  the  second 
step. 

06.  The  use  of  either  forward  or  backward  spatial  differences  to 
represent  the  advective  terns  of  the  transport  equation  is  closely  re¬ 
lated  to  the  characteristics  of  the  hyperbolic  equation.  Consider  the 
case  of  a  =  0  in  Equation  26.  The  transport  equation  then  states 
that  D$/Dt  =  0  along  the  characteristic  direction  given  by  dx/dt  =  u 
Therefore,  from  the  illustration  below 


as  shown;  whereas,  if  u  is  negative,  it  falls  between  i  and  i  +  1 
on  the  n  -  1  level.  The  major  problem  in  determining  <j>?  is  to  deter¬ 
mine  <t>#  ^  .  If  linear  interpolation  is  used  between  i  and  i  -  1  , 
the  resulting  expression  for  corresponds  to  the  use  of  a  forward- 

in-time  ar 1  oackward-in-space  representation  of  the  basic  equation. 
However,  i  u  is  negative  and  a  linear  interpolation  between  i  and 
i  +  1  is  employed,  the  resulting  expression  is  equivalent  to  the  use  of 
forward  differences  in  the  spatial  derivative.  If  a  linear  interpolation 
from  i  -  1  to  i  +  1  is  used,  centered  differences  result.  In  addi¬ 
tion,  one  could  use  higher  order  interpolating  schemes  such  as  a  qua¬ 
dratic  polynomial  to  interpolate  between  i  -  1  ,  i  ,  and  i  +  1  which 
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yields  Leith's  scheme  (Roaehe  1972).  This  scheme  is  0(At“,  )  ,  ever, 

though  only  two  time  levels  appear  in  the  difference  equation.  It  can 
be  shown  that  this  approach  is  equivalent  to  replacing  trie  second  order 
time  derivative  in  the  series  expansion  by  a  second  order  spatial  deriva¬ 
tive,  as  previously  discussed. 

07.  One-sided  differences,  i.e.,  the  forward  or  backward  schemer;, 
introduce  an  artificial  dissipation  into  the  solution  similar  to  the 
case  where  a  f  0  ;  whereas,  centered  differences  do  not  introduce  this 
dissipation.  However,  the  one-sided  differences  preserve  what  is  Labeled 
by  Roaehe  (1972)  as  the  "transportive  property,"  which  is  not  the  case 
with  centered  differences.  The  transportive  property  is  related  to 
whether  the  parameter  4>  is  numerically  aavected  solely  in  the  direc¬ 
tion  of  the  flow,  as  theoretically  it  should  be. 

88.  In  the  space  differencing  discussed  above,  only  first  or 

second  order  schemes  have  been  discussed.  However,  higher  order  spatial 

schemes  can  be  developed  and  have  been  utilized,  in  particular  in  the 

work  of  Abbott  (1979).  In  the  DHI  models  (Hinstrup  1977),  Everett's 

12-point  interpolating  polynomial  in  two  dimensions  is  used  to  generate 

a  fourth  order  transport  scheme  that  conserves  mass,  advects  correctly 

the  center  of  mass,  i.e.  ,  maintains  the  transportive  property,  has  no 

2  2 

artificial  dispersion  (proportional  to  3  <p/dx  ),  and  in  addition  con¬ 
serves  third  and  fourth  moments  of  the  distribution  of  $  .  A  disadvan¬ 
tage  of  such  higher  order  schemes  that  extend  over  several  grid  points 
is  the  difficulty  encountered  near  boundaries. 

89.  Holly  and  Preissman  (1977)  present  a  method  of  constructing 

higher  order  schemes  that  utilize  only  two  grid  points.  Their  method 
ce  ..  -ound  the  use  of  Hermitian  interpolating  polynomials  rather 

than  interpolating  polynomials  that  extend  over  several  net  points . 
Hermitian  polynomials  are  constructed  such  that  not  only  the  function 
but  also  derivatives  of  the  function  are  required  to  satisfy  known  con¬ 
ditions  at  only  two  points.  Numerical  schemes  based  on  this  concept  are 
referred  to  as  "two-point  higher  order"  methods  to  emphasize  the  fact 
that  by  using  function  derivatives,  one  can  obtain  higher  order  one- 
dimensional  schemes  using  information  at  only  two  points.  In  fact,  the 
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authors  indicate  that  results  from  comparative  studies  show  that  the  use 
of  derivatives  to  obtain  a  third  degree  interpolating  polynomial  is  more 
accurate  than  one  using  a  third  or  fourth  degree  polynomial  based  on 
additional  points. 

90.  Such  a  method,  of  course,  requires  both  the  dependent  variable 
and  its  space  derivative  as  initial  and  boundary  conditions.  However, 
through  an  example  application ,  Holly  and  Preissman  (1977)  show  that  the 
inconsistencies  introduced  between  trie  dependent  variable  and  its  deriva¬ 
tives  as  estimated  from  initial  given  values  of  the  variable  will  have  a 
minor  influence  on  the  results.  Although  an  extension  of  the  method  to 
two  dimensions  is  not  presented,  some  preliminary  computational  results 
are.  The  authors  indicate  that  such  an  extension  to  two  dimensions  pre¬ 
serves  the  favorable  accuracy  characteristics  observed  in  one  dimension. 

Consistency,  con- 
vergence,  and  stability 

91.  A  finite  difference  scheme  is  said  to  be  consistent  if  when 
one  expands  the  discrete  system  in  Taylor's  series  form  by  retaining  the 
higher  order  terms,  all  the  terms  of  the  differential  equation  (with 
possible  additional  terms)  are  generated.  In  addition,  in  the  limit  as 
the  time  and  spatial  steps  approach  zero  independently,  all  of  the  addi¬ 
tional  terms  must  go  to  zero. 

92.  In  order  for  a  numerical  solution  to  be  meaningful,  it  must 
be  a  good  approximation  of  the  exact  solution  of  the  differential  equa¬ 
tions.  Convergent  finite  difference  schemes  are  those  for  which  the 
solution  of  the  difference  equations  converges  to  the  exact  solution  as 
the  size  of  time  and  spatial  steps  approach  zero.  The  convergence  of 
finite  difference  solutions  of  the  nonlinear  equations  governing  fluid 
motions  cannot  be  proved  analytically,  and  thus,  one  must  resort  to  the 
use  of  intuition  or  preferably  a  comparison  of  numerical  results  with 
laboratory  and/or  field  data  to  demonstrate  that  the  numerical  scheme 
does  indeed  model  the  physical  processes  represented  mathematically  by 
the  governing  differential  equations. 

93.  In  a  rigorous,  mathematical  sense,  a  finite  difference  scheme 
is  stable  if  two  solutions  that  are  arbitrarily  close  to  each  other  at  a 
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given  time  remain  arbitrarily  close  for  all  time.  In  a  practical  sense, 
one  considers  a  particular  scheme  stable  if  the  solutions  do  not  grow 
unbounded.  For  economic  reasons,  in  the  numerical  calculation  of  space- 
and  time-dependent  hydrodynamic  problems,  one  desires  to  use  as  large  a 
space  and  time  step  as  possible  and  still  obtain  the  desired  level  of 
accuracy  and  physical  detail.  However,  in  addition  to  these  restric¬ 
tions,  the  stability  of  the  finite  difference  scheme  dictates  the  size 
of  the  integration  difference  steps  that  can  be  employed. 

9b.  Explicit  finite  difference  schemes  are  conditionally  stable; 
i.e.,  stable  computations  will  result  so  long  as  the  space  and  time  steps 
satisfy  what  are  known  as  "stability  criteria."  In  free  surface  hydro- 
dynamic  modeling,  the  most  severe  of  these  criteria  is  usually  the 
Courant  condition  on  a  gravity  wave. 

At  < 

^ih 

which  states  that  the  time  and  spatial  steps  are  restricted  such  that  a 
gravity  wave  will  not  propagate  over  more  than  one  spatial  step  within 
the  prescribed  time  step.  Additional  stability  criteria  presented 
below 


At  <  Ax/u 
At  <  Az2/2A 
At  <  gh 

are  related  to  the  velocity  of  a  fluid  particle,  the  rate  of  diffusion, 
and  the  speed  of  internal  waves,  respectively. 

95*  All  or  some  of  these  restrictions  may  be  eliminated  by  various 
finite  difference  schemes.  For  example,  fully  implicit  schemes  can  be 
constructed  that  are  unconditionally  stable,  at  least  in  a  linear  sense; 
whereas,  mixed  implicit-explicit  schemes,  such  as  that  of  Edinger  and 


Buenak  (1979)  ,  may  be  constructed  to  remove  one  or  more  of  the  more 
severe  criteria  while  continuing  to  be  restricted  by  the  less  severe 
ones.  Each  finite  difference  scheme  has  its  own  advantages  and  diffi¬ 
culties,  and  which  scheme  is  best  often  depends  upon  the  particular 
problem.  For  example,  one  may  be  able,  from  a  stability  standpoint,  to 
use  an  unlimited  time  step  in  an  implicit  scheme  as  opposed  to  perhaps 
a  rather  small  time  step  in  an  explicit  scheme.  However,  if  the  physi¬ 
cal  character  of  the  problem,  such  as  rapidly  varying  input  boundary 
conditions,  forces  tiie  use  of  a  relatively  small  time  integration  step 
in  the  implicit  code,  one  may  find  that  ar.  explicit  model  is  actually 
more  economical  due  to  the  simplicity  of  the  solution  technique. 

96.  Stability  of  a  finite  difference  scheme  car  be  related  to  the 
concept  of  artificial  viscosity  or  diffusivity,  whioa  has  been  previously 
discussed.  Using  Hirt's  method  of  analysis,  as  opposed  to  th=  more  elab¬ 
orate  von  Neumann  analysis  in  which  the  growth  of  a  Fourier  component  is 
investigated  (see  Roache  1972),  consider  the  stability  of  a  forward-in¬ 
time  and  centered-in-space  representation  of  Equation  2 6: 
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Replacing  the  discrete  values  above  by  a  Taylor  series  expansion  and 
making  use  of  the  initial  differential  equation  yields 
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(39) 


It  can  be  seen  that  as  At  and  Ax-K)  ,  the  above  equation  reduces  to 
the  original  differential  equation;  therefore,  the  difference  scheme  is 

consistent.  However,  it  will  not  be  stable  unless  the  effective  dissi- 

2 

pative  coefficient  a  -  u  At/2  is  greater  than  zero,  since  the  physical 
nature  of  such  a  coefficient  is  to  smear  a  disturbance.  Thus,  a  negative 
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coefficient  is  a  physical  impossibility.  The  term  -u  At/ 2  is  referred 
to  as  the  "negative  artificial  viscosity"  or  "diffusitivity"  of  the 
scheme.  One  can  now  see  why  a  forward  time  integration  with  centered 
spatial  derivatives  will  result  in  a  completely  unstable  scheme  when 
applied  to  pure  hyperbolic  equations,  i.e.,  u  =  0  in  Equation  2 6. 

If  backward  differences  (u  >  0)  are  used  to  replace  the  spatial  deriva¬ 
tive,  the  effective  dissipative  coefficient  becomes 

a  =  a  +  uAx/2  -  u^At/2 
e 

Therefore,  one-sided  differences  for  spatial  derivatives  increase  the 
effective  dissipative  coefficient,  which,  results  in  a  more  stable, 
although  theoretically  less  accurate,  scheme. 

* 


\ 
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I  ALT  IV:  THREE-DIMENSIONAL  HYDRODYNAMIC  MODELS 

97-  A  relatively  wide  range  of  numerical  three-dimensional  hydro- 
dynamic  models  currently  exist.  All  that  have  been  investigated  utilize 
finite  differences  to  obtain  numerical  solutions  of  the  governing  equa¬ 
tions  and  furthermore  all  employ  explicit  time  differencing.  The  gen¬ 
eral  opinion  in  the  past  concerning  the  use  of  3-D  implicit  schemes  has 
been  that  due  to  the  extremely  large  matrices  that  have  to  be  inverted 
for  a  completely  implicit  model,  such  schemes  would  require  excessive 
computing  time.  In  addition,  apparently  schemes  such  as  Leendertse 
(1967)  employs  in  his  two-dimensional  work  (alternating  direction 
implicit — ADI)  have  not  been  used  for  various  reasons.  First,  such 
schemes  require  all  computational  arrays  to  be  in  the  computer's  fast 
memory  for  efficient  computation,  which  would  put  a  considerable  re¬ 
straint  on  the  array  sizes  of  a  three-dimensional  model.  In  addition, 
such  schemes  place  restrictions  on  the  formulation  of  the  finite  dif¬ 
ference  representation  of  various  terms  in  the  equations. 

98.  The  time-dependent  and  variable  density  three-dimensional 
models  of  Simons  (1973),  Lick  (1976),  Leendertse  et  al.  (1973),  Waldrop 
and  Tatom  (1976),  and  Spraggs  and  Street  (1975)  are  discussed  in  some 
detail  below.  Other  less  general  three-dimensional  numerical  hydro- 
dynamic  models  exist,  such  as  those  of  Gedney  and  Lick  (1970) »  Liggett 
(1970),  and  Bonham-Carter  et  al.  (1973).  However,  for  the  computation 
of  flows  in  stratified  reservoirs,  only  those  models  that  are  time- 
dependent  and  allow  for  a  variable  density  are  of  interest. 

Simons'  3-D  Lake  Model 

99.  The  modeling  of  stratified  fluid  flow  may  be  accomplished  in 
two  ways:  (a)  a  layered  model  in  which  the  fluid  is  made  up  of  discon¬ 
tinuous  layers  within  which  all  fluid  properties  such  as  density  and 
viscosity  are  uniform  and  (b)  a  continuous  model  in  which  the  density 
is  varied  continuously.  Historically,  in  numerical  models  developed 

by  meteorologists  and  oceanographers,  the  three-dimensional  model  has 


57 


been  viewed  as  a  superposition  of  layers  of  fluid  separated  by  material 
interfaces.  The  reason  for  this  is  partly  physical,  since  during  cer¬ 
tain  periods  a  body  of  water  may  become  so  stratified  that  strong,  den¬ 
sity  discontinuities  can  exist.  Jn  the  other  hand,  if  the  vertical 
resolution  of  the  model  is  sufficiently  large,  any  type  of  stratifi¬ 
cation  can  be  handled  by  a  straightforward  three-dimensional  finite 
difference  grid,  i.e.,  a  sequence  of  rigid  permeable  horizontal  levels. 

100.  Simons'  (1973)  model  is  a  multilayered  model,  which  employs 
the  principles  and  terminology  of  layered  models  while  retaining  the 
capability  of  treating  the  layers  as  being  separated  by  permeable 
rigid  interfaces  (either  horizontal  or  sloping)  as  well  as  treating  the 
interfaces  in  the  usual  layered  manner  as  moving  material  surfaces.  The 
equations  for  the  layered  system  are  obtained  by  vertical  integration  of 
the  governing  equations  (written  in  the  conservative  form)  over  each 
layer  as  opposed  to  applying  the  equations  at  given  levels  and  replacing 
the  vertical  derivatives  by  finite  differences.  The  primary  dependent 
variables  are  the  layer  thickness  or  vertical  velocity  and  the  layer- 
averaged  horizontal  velocity  components  as  well  as  the  temperature. 

101.  Simons  invokes  the  Boussinesq  approximation  and  assumes  that 
vertical  accelerations  are  negligible;  i.e.,  the  pressure  is  hydrosta¬ 
tic.  With  the  assumption  of  the  Boussinesq  approximation,  the  equation 
of  mass  continuity  reduces  to  the  incompressibility  condition,  which 
implies  that  the  vertical  fluid  motion  is  directly  related  to  the  diver¬ 
gence  of  the  horizontal  flow.  With  the  hydrostatic  pressure  assumption 
replacing  the  vertical  momentum  equation,  the  vertical  component  of 
velocity  cannot  be  computed  in  the  s;nne  manner  as  the  two  horizontal 
components.  Instead,  the  equation  expressing  incompressibility  is 
integrated  over  a  layer  to  yield  an  equation  whose  primary  purpose  is 

to  compute  water  displacements  from  a  given  distribution  of  horizontal 
velocities.  From  this  equation,  one  can  determine  either  the  displace¬ 
ment  of  a  material  surface  or  the  vertical  velocity  of  the  fluid  through 
a  rigid  interface,  if  given  the  appropriate  boundary  conditions  at  the 
free  surface,  at  the  interface,  and  at  the  bottom.  The  computation 
starts  with  the  bottom  layer  and  proceeds  upward. 


102.  The  equation  of  state  is  such  that  the  density  is  assumed 
to  be  a  quadratic  function  of  the  temperature.  An  eddy  coefficient 
model  is  employed  to  approximate  the  exchange  of  energy  between  the 
large-scale  flow  and  the  smaller  turbulent  eddies .  It  appears  that 
constant,  but  different,  horizontal  and  vertical  coefficients  are 
employed.  Simons  does  indicate  that  the  vertical  eddy  diffusivity 
depends  on  the  static  stability  3p/3z  of  the  water  column,  and  he 
allows  it  to  attain  very  large  values  for  unstable  situations  in  order 
to  simulate  the  net  effects  of  convective  overturning. 

103.  The  time  integration  scheme  employed  by  Simons  uses  centered 
differences,  where  the  pressure  gradient  terms,  the  divergence  terms, 
the  Coriolis  terms,  and  the  nonlinear  terms  are  evaluated  at  a  time  step 
centered  between  the  old  and  new  time,  while  the  dissipative  and  diffu¬ 
sion  terms  are  evaluated  at  the  old  time  step.  Centered  differences  are 
also  used  to  replace  spatial  derivatives,  and  thus,  the  finite  differ- 
ence  model  is  almost  0(At  ,  Ax  )  .  Linear  interpolation  is  used  to 
provide  values  of  variables  at  points  where  they  are  not  defined. 

IOU.  An  interesting  aspect  of  Simons'  model  is  his  use  of  two 
different  time  steps.  The  surface  and  internal  computations  are  de¬ 
coupled  such  that  a  small  time  step  governed  by  the  Courant  condition 
is  used  to  compute  the  water  surface  elevations;  whereas,  a  much  larger 
time  step  governed  by  the  speed  of  a  fluid  particle  is  used  to  compute 
the  internal  flow  and  the  temperature.  This  is  accomplished  as  follows. 
The  layer-averaged  equations  are  added  to  create  a  vertically  averaged, 
i.e.,  one-layer,  model  for  gross  fluid  flows,  which  are  then  used  to 
drive  the  free  surface.  The  layer-averaged  equations  then  use  the  re¬ 
sults  of  the  vertically  averaged  model  to  produce  the  internal  flow 
field.  Simons  indicated  that  in  Lake  Ontario,  with  a  grid  mesh  of  5  km, 
the  surface  elevation  and  the  vertically  integrated  flow  were  computed 
with  a  time  step  of  the  order  of  one  minute,  while  the  internal  flow  and 
temperature  were  predicted  with  a  time  step  of  the  order  of  30  min. 

Leendertse's  3-D  Estuary  Model 

105.  This  variable  density  model  (Leendertse  1973)  has  been 
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developed  for  the  computation  of  the  hydrodynamics  of  estuaries  and 
coastal  seas.  Some  of  the  assumptions  and  formulations  are  not  direct! 
applicable  to  freshwater  reservoirs  where  density  effects  are  due  to 
temperature  variations  rather  than  salinity;  however,  it  is  believed 
that  such  an  extension  would  not  be  difficult. 

106.  The  approach  taken  in  the  formulation  of  tiie  equations  is 
similar  to  Simons  (1973)  in  that  tne  basic  three-dimensional  equations 
are  integrated  over  a  vertical  layer  to  yield  layer-averaged  equations. 
Unlike  Simons,  however,  Leendertse's  model  does  not  allow  movable  mate¬ 
rial  interfaces.  The  water  body  is  represented  by  rigid  permeable  hori 
zontal  surfaces  with  the  thickness  of  each  interior  layer  constant  in 
space  and  time,  although  the  thickness  of  each  layer  is  not  necessarily 
the  same.  The  top  layer  that  contains  the  surface  is,  of  course,  rep¬ 
resented  by  a  time-varying  and  spatially  varying  thickness. 

107.  The  basic  three-dimensional  equations  are  written  in  the 
conservative  form  to  insure  that  mass,  momentum,  etc.,  are  neither 
created  nor  destroyed  by  the  computational  scheme.  Before  the  layer 
integration  is  performed,  the  Boussinesq  approximation  is  assumed  and 
the  pressure  is  assumed  to  be  hydrostatic.  Therefore,  as  in  Simons' 
model,  the  vertical  component  of  the  fluid  velocity  must,  be  computed 
from  the  layer-averaged  condition  of  incompressibility. 

108.  Approximate  eddy  viscosity  models  that  consider  only  the 
diagonal  components  of  the  viscosity  tensor  are  employed  to  represent 
the  subgrid-scale  motions.  The  momentum  and  mass  dispersion  coeffi¬ 
cients  are  assumed  to  be  constant  in  the  horizontal  dimensions  of  the 
flow,  although  they  can  differ  in  the  two  directions..  The  vertical 
exchange  coefficients  are  calculated  with  a  mere  sophisticated  model 
that  takes  into  account  the  vertical  velocity,  the  concentration  tra- 
dient,  and  the  stability  of  the  flow  according  to  the  hi ohardson  number 

109.  Since  the  model  was  developed  for  an  estuarine  or  const  al 
environment,  the  equation  of  state  relates  the  fluid,  density  to  the 
salinity.  If  the  model  were  to  be  applied  to  a  fir:;;ivi'sr  env .  r  ...ament , 
a  new  equation  of  state  relating  the  density  to  temperature  won!  i,  of 
course,  need  to  be  substituted.  Tn  addition,  surface  heat  exchnn  -<■ 


would  have  to  be  accounted  for  through  the  surface  boundary  condition. 

110.  At  the  boundaries  of  the  water  body  to  be  computed,  all 
diffusion  coefficients  are  set  to  zero,  as  are  the  velocities  perpen¬ 
dicular  to  the  boundary.  In  this  manner,  no  mass  fluxes  or  diffusive 
transports  of  salt  result.  At  the  surface,  the  boundary  stress  due  to 
the  wind  is  computed  from  a  quadratic  law.  A  similar  quadratic  expres¬ 
sion  employing  the  Chezy  coefficient  is  used  to  represent  the  dissipa¬ 
tion  of  momentum  at  the  bottom  through  the  bottom  shear  stress.  When 
employing  the  layer-averaged  approach,  interfacial  shear  stress  terms 
show  up  in  the  resulting  equations  for  the  layer-averaged  variables. 

As  for  the  boundary  stress  specifications,  a  quadratic  relationship 
between  the  interlayer  stresses  and  the  velocity  differences  of  adjacent 
layers  is  assumed  applicable. 

111.  The  spatial  grid  used  in  the  finite  difference  formulation 
is  similar  to  that  employed  by  Leendertse  (1967)  in  his  two-dimensional 
work  where  velocities  are  defined  on  the  faces  of  a  cell.  However,  the 
water  surface  elevations,  which  are  determined  from  an  equation  obtained 
by  summing  the  layer-averaged  incompressibility  equation  over  the  water 
column,  are  defined  at  the  corners  of  the  top  layer  of  cells  rather  than 
at  the  cell  center  as  in  the  2-D  model.  Pressure,  density,  and  salinity 
are  defined  at  cell  centers. 

112.  Centered  differences  are  used  for  both  time  and  spatial 

integrations.  Therefore,  the  resulting  finite  difference  scheme  is 
2  2 

0(At  ,  Ax  )  except  that  the  diffusion  terms  are  taken  at  the  lower 
time  level,  i.e.,  t  =  (n-l)At  ,  since  otherwise  the  computation  becomes 
unstable.  Since  centered  differences  are  used  to  replace  the  time  de¬ 
rivatives,  initial  information  at  two  time  levels  is  required.  To  re¬ 
duce  the  time-splitting  tendency  of  such  a  scheme,  a  single  forward 
differencing  step  is  used  to  obtain  initial  information  on  the  second 
time  step. 

113.  Since  an  explicit  time  integration  scheme  has  been  utilized, 
the  basic  stability  criterion  is  once  again  the  Courant  condition. 
Although  the  computations  with  the  adopted  scheme  are  extensive, 
Leendertse  indicates  that  the  model  is  well  within  the  range  of 
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i  i  ?k*  s  Thermal  i  lame  Molel 


•  I...  a  variabi ••  iensity,  heat-conducting  model  for  studying  the 
field  nurr.:.;::  ii.v  r.ue  :•  unt  if  iljoha r.-;e  of  a  river  or  power  plant 
int  •  a  tody  of  water  -.a.  teen  ievelopel  by  Lick  U  T.C)  at  lane  Western 


11'  .  Panic  an.;  tior.n  made  :  n  ihe  model  development  are  that  tne 
boussinenq  approx imai. i cm  holds  and  that  tne  pressure  is  nydrostat ic.  «c 
previously  dir.can.ie  i,  the  hy  irontatic  pressure  assur.pt ion  reduces  the 
order  of  the  system  of  equations  and  the  computational  effort  required. 

A  major  difference  between  the  Lick  model  and  other  t. hre diraens i ona* 
models  studied  is  the  assumption  of  a  rigid  lid  at  the  surface.  Inis 
approximation,  which  removes  tiie  surface  gradient  terms  from  the  momen¬ 
tum  equations,  is  made  to  eliminate  surface  gravity  waves  and  tne  smu: 
time  scales  associated  with  them.  Inis,  of  course,  greatly  increases 
the  maximum  allowable  time  step  without  encountering  stability  ; rot lem.  . 
Witt;  the  rigid-lid  approximation,  the  limiting  time  step  in  no  longer 
given  by  the  extremely  restrictive  gravity  wave  Courant  coalition  big. 
instead  is  usually  fixed  by  the  speed  of  a  water  particle,  i.e.. 

At  <  Ax/u 

or  perhaps  by  the  speed  of  an  internal  wave,  i.e.. 


At  < 


Ax 


rigid-lid  assumption,  the  pressure  is  no  longer 


•it:-.  :.(>!  Le  m*  *. .surface  bus  instead  varies  with  the  internal  flow 

fit  L;.  ,  a:,  additional  equation  must  be  derived  and  subsequently 

solve  1  ft!  yieli  tlae  .surface  pressure.  Once  the  surface  pressure  is 
ie*-  iutmul  pressure.-  are  determined  from  the  hydrostatic 

1*7.  The  :  Is  son  e  puat ion  for  the  surface  pressure  is  derived  by 
taking  t.ne  divergence  of  tne  vertically  integrated  momentum  equations 
and  making  use  of  tne  vertically  integrated  continuity  and  ’hydrostatic 
pressure  equations.  The  fora  of  this  equation  becomes 


F  ( u ,  v ,  w ,  T ) 


where  r  is  the  surface  pressure.  Lick  indicates  that  even  though  a 
rigid  lid  has  been  assumed,  one  can  compute  surface  displacements  (ne¬ 
glect  ins  the  transient  motion  due  to  gravity  waves)  by  interpreting  the 
surface  pressure  as  a  pressure  due  to  a  height  of  water  above  or  below 
the  location  of  the  rigid  lid.  The  numerical  solution  of  the  pressure 
Poisson  equation  at  each  time  step  is  accomplished  by  using  the  ADI 
method . 

118.  The  diagonal  components  of  the  eddy  coefficient  tensors  are 
used  to  account  for  the  turbulent  subgrid-scale  motions.  The  horizontal 
eddy  coefficients  are  assumed  constant,  but  the  vertical  eddy  coeffi¬ 
cient  are  a  function  of  the  temperature  gradient  and  other  parameters. 
This  lependenee  on  temperature  is  given  by 


A 

v 


a 


where  A  is  the  vertical  eddy  di ffusivity  and  a  and  6  are  con- 
v 

stunt  depending  on  the  local  conditions.  The  constant  a  is  equal  to 
the  vertical  eddy  diffusivity  under  neutral  stability  conditions.  Typi¬ 
cal  values  for  a  and  6  are  70  cm" /sec  and  200  cm  /°C-sec,  respec¬ 
tively.  Lick  handles  a  static  instability  in  the  same  manner  as  Pinions; 
i.e.,  extensive  mixing  is  assumed. 
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■j  119*  An  interesting  aspect  of  Lick's  model  is  his  use  of  a  linear 

|  transformation  in  the  vertical  direction  such  that  the  z  ,  y  ,  z 

Cartesian  system  is  transformed  to  a  x  ,  y  ,  a  system  where  1  <.  7 
'  <_  1  ,  with  the  bottom  corresponding  to  a  =  0  and  the  top  to  e  =  1  . 

i  The  o  coordinate  is  defined  by 


°  =  1  +  h 


where  h(x,  y)  is  the  depth  of  the  water  body.  With  sue:;  a  trans 
tion,  irregular  bottoms  can  be  handled  more  accurately  and  efficit 
For  example,  the  usual  finite  difference  representation  of  .u 
bottom  of  the  IRK  flume  (applications  to  be  discussed  later)  is  a. 
below. 


£3^ 


W‘ 


However,  the  use  of  the  transformation  above  vields  :..s-  following 
physical  representation,  with  the  actual  numerical  computations  being 
performed  on  a  rectangular  transformed  grid: 


:*  sr.oul  i  be  noted,  however,  taut  since  the  t r an s format  Ion 
>r::;ai ,  wh>*n  transf  UT.ir.g  the  governing  equations,  t  •••i.n 
second  iffivat ;  vo  diffusive  terms  contain  cross-  :•••*!  vu l  ivs 
s:  at  ini  'O' ir  1  i  nates. .  Lick  a.  .  one.-  goat  ti.o  1 1  f  fus  i  v<  • 
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derivatives  of  the  depth  are  negligible  with  respect  to  those  diffusive 
terns  containing  only  the  depth  and  thus  drops  such  ferns. 

120.  Solid  boundaries  are  taken  as  no-slip,  impermeable ,  adia¬ 
batic  surfaces.  A  heat  transfer  condition  proportional  to  the  temper¬ 
ature  difference  (surface  temperature  minus  equilibrium  temperature) 
and  a  wind-dependent  stress  are  imposed  at  the  surface.  "ormal  deriva¬ 
tive  pressure  boundary  conditions  are  derived  from  the  appropriate 
vertically  integrated  momentum  equation.  At  open  boundaries  far  from 
the  point  of  discharge  of  the  river  or  plume,  the  normal  derivatives 

of  the  velocities  and  the  temperature  are  zero. 

121.  As  opposed  to  the  layer-averaged  approach  of  Simons  and 

Leendertse,  Lick  performs  a  straightforward  finite  differencing  of  the 
governing  3-D  equations.  A  forward  time-differencing  integration  scheme 
is  utilized  along  with  centered  differences  for  the  spatial  derivatives. 
Thus,  the  finite  difference  scheme  is  of  the  first  order  in  time  and 
second  order  in  space,  i.e.,  0(At,  Ax  )  .  The  computational  grid  is 

such  that  the  horizontal  velocity  components,  u  and  v  ,  are  defined 
as  the  cell  corners  with  the  vertical  component  defined  at  the  middle 

of  the  top  and  bottom  face  of  the  cell.  The  pressure  and  temperature 
are  defined  at  the  cell  center,  except  for  the  surface  pressure,  which 
is  computed  at  the  center  of  the  top  face. 


Waldrop-Tatom  3-D  Plume  Model 

122.  There  are  actually  two  versions  of  this  extremely  versatile 
three-dimensional  variable  density  model  (Waldrop  and  Tatom  1976).  One 
employs  she  hydrostatic  pressure  assumption,  and  the  other  retains  the 
complete  vertical  momentum  equation.  Both  utilize  the  Boussinesq  ap¬ 
proximation  and  both  neglect  Coriolis  effects.  It  appears  from  Waldrop 
and  Tatom  (1976)  that  the  hydrostatic  pressure  version  solves  the  non¬ 
conservative  form  of  the  basic  governing  equations;  whereas,  Tatom  and 
Smith  (1979a)  indicate  that  the  conservative  form  of  the  equations  are 
solved  in  the  version  that  does  not  make  the  hydrostatic  assumption. 
Both  versions  solve  the  governing  equations  transformed  into  orthogonal 

6  5 


curvilinear  coordinates.  This,  of  course,  allows  for  store  accurate 
inodelin;'  of  curved  boundaries  such  as  river  bends. 

If 3.  In  the  nonhydrostatic  version,  the  pressure  is,  written  nr, 
the  sura  of  the  hydrostatic  pressure  and  the  dynamic  pressure  1/1  f  v'  , 
and  a  Poisson  equation  for  the  dynamic  pressure  is  derived  for  solution 
over  the  complete  3-1'  field.  The  Richardson  iterative  technique  is  em¬ 
ployed.  It  raight  be  noted  that  in  the  Poisson  pressure  equation,  terms 
that  involve  the  horizontal  density  gradients  have  been  negloetei.  how¬ 
ever,  it  does  appear  that  horizontal  density  gradients,  as  reflected 
through  the  hydrostatic  component  of  the  pressure,  are  included  in  the 
velocity  computations  from  the  momentum  equation. 

12^.  With  the  retention  of  the  complete  vertical  momentum  equa¬ 
tion,  a  fully  convective  model  that  can  handle  buoyancy  effects,  i.e., 
unstable  density  profiles,  is  realized.  The  vertical  component  of  the 
velocity  is  now  determined  from  the  vertical  momentum  equation  as 
opposed  to  its  solution  from  the  incompressibility  condition  in  the 
hydrostatic  version. 

125.  The  effects  of  turbulence  are  included  through  the  use  of 

eddy  coefficients.  The  horizontal  eddy  viscosity  coefficient  £„  is 

n 

derived  from  a  mixing  length  equation  for  open  channels  in  the  form 
eH  =  0.l6(z  -  zB)2  U  -  z)/U  -  z3)  1 3-^u2  +  v^/Szj 


where 

z  =  z_  at  the  bottom 

D 

z  =  £  at  the  free  surface 

which  provides  the  largest  values  of  e  in  deep  regions  with  large 
velocity  gradients  in  the  vertical;  whereas,  in  shallow  and/or  low  flow 
regions  e  is  small.  The  horizontal  eddy  diffusivity  Au  is  related 

n  n 

to  the  eddy  viscosity  by 


Ag  =  1.33eh 


126.  In  turbulent  flows,  density  stratification  inhibits  the 


vertical  exchange  of  both  heat  and  momentum  as  well  as  the  mass  of  any 
constituent.  The  V.'aldrop-Tatoin  models  allow  the  vertical  eddy  coeffi¬ 
cients  to  be  functions  of  the  stratification  through  their  dependence- 
on  the  Richardson  number  in  the  following  manner: 

-1  /? 

e  =  e..  (1  +  10R.  )  ' 

v  H  1 


Av  =  (1  +  3.33PL) 


-3/2 


where 


fh  =  -(g/p)(3p/3z) 


(a  \Ju2  +  v2/3z) 


It  might  be  noted  that  although  the  eddy  coefficients  are  allowed  to 
vary  spatially,  spatial  derivatives  of  the  coefficients  have  been 
neglected  in  the  model. 

127.  At  solid  boundaries,  reflection  boundary  conditions  are 
imposed  to  simulate  slip  boundaries.  Therefore ,  with  solid  walls  as¬ 
sumed  to  lie  between  the  last  two  grid  points,  fictitious  values  of 
dependent  variables  on  the  opposite  side  of  a  wall  are  set  to  prevent 
mass,  momentum,  or  energy  transfer  through  the  boundaries.  Velocities 
normal  to  the  wall  are  set  as  the  negative  of  the  value  immediately 
inside  in  order  to  make  the  normal  velocity  zero  at  the  wall,  but  the 
tangential  component  is  set  equal  to  its  value  inside  since  with  slip 
walls,  the  wall  does  not  influence  the  tangential  flow.  Derivatives 

of  the  temperature  normal  to  solid  walls  are  set  equal  to  zero  to  insure 
no  transfer  of  heat. 

128.  The  velocity  profile  near  the  bottom  is  assumed  to  be  loga¬ 
rithmic.  Thus,  the  equation  below  is  used  to  help  set  the  horizontal 
velocity  components  at  all  grid  points  adjacent  to  the  bottom  in  the 
solution  of  the  momentum  equations: 


1 


u  =y{TJ-p 


z 


k 


+  8.5 


where 


T  =  snear  stress 
o 

(z  -  z_)  =  height  above  bottom 

D 

k  =  diameter  of  the  average  roughness 
The  actual  values  chosen  are  such  that  the  finite  difference  representa¬ 
tion  of  the  velocity  gradients  9u/9z  and  9v/9z  near  the  bottom  match 
the  gradient  specified  by  the  equation  above-  As  noted,  this  is  the 
procedure  for  determining  bottom  velocities  for  use  in  the  momentum 
equations.  However,  in  the  transport  equation  for  temperature,  the 
velocities  at  points  adjacent  to  the  bottom  are  determined  from  an 
actual  fit  of  the  logarithmic  profile  rather  than  by  forcing  the  proper 
gradient.  In  the  computation  of  the  free  surface,  a  control  volume  is 
formed  between  the  top  grid  plane  and  the  free  surface.  Since  the  three 
velocity  components  from  previous  computations  at  a  particular  time  line 
are  known,  the  mass  transported  into  and  out  of  the  control  volume  can 
be  computed.  The  free  surface  is  then  adjusted  to  insure  conservation 
of  mass.  In  the  current  versions  of  the  model,  the  time  integration  is 
essentially  a  forward  difference,  but  with  an  additional  step  that 
Waldrop  and  Tatom  (1976)  indicate  helps  to  stabilize  the  computations. 
This  is  accomplished  with  the  following  scheme: 


where  (9u/9t)n  ^  is  saved  from  computations  at  the  previous  time  step. 
It  would  appear  that  this  is  equivalent  to  replacing  the  time  derivative 
at  t  =  (n  -  l/2)At  by  a  forward  difference  between  (n  +  l)At  and 
nAt  .  Thus,  the  scheme  is  still  only  first  order  in  time.  Centered 
differences  are  used  in  the  diffusive  terms,  while  one-sided  windward 
(either  forward  or  backward,  depending  upon  the  direction  of  flow) 
differences  are  used  in  the  finite  difference  representation  of  the 
advective  terms.  Thus,  it  would  appear  that  the  solution  scheme  is 
0(At,  Ax)  . 

129*  As  noted  in  a  previous  section,  the  computational  grid  is 
such  that  all  variables  are  defined  at  the  same  point,  with  uneven 
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spacing  of  those  points  allowed  for  more  flexible  resolution.  Waldrop 
and  Tatom  indicate  that  a  transformation  of  the  x  ,  y  ,  z  coordinates 
such  that  even  increments  in  the  transformed  system  produce  uneven  spac¬ 
ing  of  the  grid  points  in  the  physical  system  is  employed.  However, 
details  of  the  transformation  are  not  discussed. 

130.  The  Waldrop-Tatom  model  is  capable  of  handling  branching 
systems  through  its  modular  concept  in  which  the  equations  are  solved 
simultaneously  in  different  branches  or  regions.  The  regions  are  con¬ 
nected  such  that  when  there  is  free  flow  between  regions,  each  region 
uses  previously  computed  information  from  the  adjacent  region  as  a 
boundary  condition.  Of  course,  the  fact  that  an  explicit  time- 
integration  scheme  has  been  employed  greatly  decreases  the  difficulty 
in  incorporating  such  a  concept.  The  handling  of  connecting  branches, 
i.e.,  connecting  regions,  in  an  implicit  model  would  be  much  more  dif¬ 
ficult  to  accomplish.  The  capability  of  handling  connecting  regions, 
allowing  for  a  variable  grid,  and  the  use  of  curvilinear  coordinates 
makes  the  Waldrop-Tatom  model  extremely  versatile. 


131.  The  nonhydrostatic  version  of  the  Waldrop-Tatom  model  and 
the  three-dimensional  model  developed  by  Spraggs  and  Street  (1975)  are 
the  only  3-D  numerical  models  studied  that  are  fully  convective  models. 
In  other  words,  the  complete  vertical  momentum  equation  is  retained  so 
that  buoyancy  effects  are  modeled  directly.  As  indicated  by  Spraggs  and 
Street,  the  primary  use  of  the  model  is  to  simulate  flows  in  which  the 
stratification  induced  by  heated  effluents  sets  up  in  a  matter  of  hours . 
No  claim  is  made  as  to  the  usefulness  of  the  present  form  of  the  model 
for  simulating  flows  over  periods  extending  over  the  time  required  for 
the  formation  of  a  natural  thermocline.  This  is  because  of  the  exces¬ 
sive  computing  time  required  due  to  the  extremely  small  time  step  im¬ 
posed  by  the  explicit  nature  of  the  solution. 

132.  As  in  the  vast  majority  of  hydrodynamic  models,  the 
Boussinesq  approximation  is  made,  which  reduces  the  conservation  of 


mass  equation  to  the  incompressibility  condition.  In  addition,  an  eddy 
viscosity  model  is  used  to  simulate  the  transfer  of  energy  from  the 
developing  flow  to  small-scale  turbulent  eddies,  i.e.,  the  subgrid-scale 
motions.  These  appear  to  be  the  only  assumptions  made  to  the  basic 
equations.  It  should  be  noted,  however,  that  one  important  restriction 
exists  in  the  basic  mathematical  development  of  the  model  due  to  the 
manner  in  which  pressure  gradients  are  handled  in  the  horizontal  momen¬ 
tum  equations. 

133.  A  reduced  pressure  P  ,  which  is  a  measure  of  the  perturba- 

K 

tions  in  the  system,  e.g. ,  caused  by  stratification  and/or  vertical  ac¬ 
celerations,  is  defined  as 

P  |p  -  V 

R'  Pr 

where  the  hydrostatic  pressure  P  is 

P,  =  ( L  -  ^  -  z )  p  g 
h  z  r 

Pr  is  the  density  of  a  reference  state,  is  a  reference  depth, 

£(x,y)  is  the  water  surface  elevation,  and  z  is  the  distance  above 
the  reference  bottom.  With  the  hydrostatic  pressure  defined  as  above 
in  terms  of  a  reference  density  that  is  not  a  function  of  (x,y),  the 
pressure  gradient  becomes 

1_  3P _ __R_  3; 

p  3x.  3x.  ®  9x. 

r  i  l 

which  does  not  allow  for  the  effect  on  the  flow  of  horizontal  gradients 
in  the  density.  It  appears  that  this  restriction  could  be  removed  by 
defining  the  hydrostatic  pressure  in  terras  of  the  spatially  varying 
density  rather  than  of  a  constant  reference  density.  Both  Edinger  and 
Buchak  (1979)  in  the  modeling  of  stratified  reservoirs  and  Hamilton 
(1975)  in  the  modeling  of  salinity-stratified  estuaries  have  indicated 
that  the  horizontal  density  gradients  are  quite  important  in  modeling, 
variable  density  flows. 
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13*+.  As  is  usually  the  case  when  the  hydrostatic  pressure  assump¬ 
tion  is  not  made,  a  Poisson  equation  for  the  reduced  pressure  P  is 

n 

derived  by  taking  the  divergence  of  the  vector  momentum  equation  and 

combining  with  the  time  derivative  of  the  incompressibility  condition. 

Derivative  boundary  conditions  on  the  pressure  at  solid  walls  are 

derived  from  the  momentum  equations;  whereas,  the  pressure  itself  is 

prescribed  at  the  free  surface.  The  solution  of  the  pressure  from  the 

three-dimensional  Poisson  equation  is  obtained  through  the  iterative 

method  called  point  Successive-Over-Relaxation  (SOR).  Spraggs  and 

Street  indicate  that  the  pressure  solution  usually  converges  within 

50  iterations.  Such  a  solution  of  a  3-D  Poisson  equation  at  each  time 

step  constitutes  a  major  portion  of  the  total  computation  time  of  the 

model.  Thus,  one  can  see  why  the  hydrostatic  pressure  assumption  has 

been  so  popular  in  the  past  in  the  development  of  hydrodynamic  models. 

135-  The  mathematical  model  is  rendered  dimensionless  through 

the  introduction  of  three  length  scales,  L  ,  L  ,  and  L  ,  such 

x  y  s 

that  any  physical  problem  is  mapped  to  the  interior  of  a  unit  cube. 

Thus,  in  the  numerical  model,  there  are  six  length  parameters — , 

Ly  ,  ;,z  ,  Ax  ,  Ay  ,  Az  .  The  first  three  are  defined  as  above, 

while  the  second  three  are  determined  by  the  number  of  computational 

cells  within  the  unit  cube.  If  L  Ax  =  L  Ay  =  L  Az  ,  the  numerical 

x  y  z 

model  is  undistorted,  and  the  computational  cells  in  the  physical 

problem  are  cubes.  Generally,  the  horizontal  length  scales  will  be 

much  larger  than  the  vertical  length  scale  giving  rise  to  a  distorted 

model  in  which  L  Ax  4  L  Ay  4-  L  Az  . 

x  y  z 

136.  The  free  surface  elevation  c,  is  computed  from  the  kine¬ 
matic  boundary  condition 


K 

3t 


w 


u 


v  IS. 

3x  dy 


where  the  vertical  coordinate  is  positive  downward.  The  solution  of 
the  free  surface  is  obtained  through  the  following  ADI  scheme,  which 
is  one  iteration  of  the  Peaceman-Rachford  scheme  with  an  acceleration 
parameter  of  1.0: 
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From  the  above  solution  technique,  it  can  be  seen  that  since  velocities 
at  the  n+1  time  level  are  required,  they  are  computed  before  the 
computations  for  the  free  surface  are  made. 

137.  The  Spraggs  and  Street  model  is  the  only  3-D  model  investi¬ 
gated  that  allows  for  tensor  eddy  coefficients,  i.e.  ,  the  off-diagonal 
terms  are  not  neglected.  The  form  of  the  eddy  viscosity  tensor  selected 

by  Spraggs  is  a  function  of  the  rate  of  strain  S  ,  i.e., 

J  mn 


e . 


ij 


=  SlAx.Ax,(S  S  ) 
1  j  mn  mn 


1/2 


where  the  Reynolds  stress  is 


u!  u'.  =  -e.  .  S.  .  (no  summation  over  i  ) 

1  J  ij  ij 

and  the  rate  of  strain  tensor  S.  .  is 

ij 

3u.  9u, 

S  =  — i.  +  — 1 
ij  3Xj  3xi 

As  Spraggs  and  Street  note,  there  is  some  question  as  to  the  value  of 
the  scaling  parameter  Q  ,  since  the  range  of  problems  that  might  be 
simulated  could  extend  from  laboratory  flume  dimensions  to  several 
hundred  kilometres  in  the  field.  A  value  of  2  =  0.01  war.  used  by 


Spraggs  and  Street  in  the  initial  testing  of  the  model.  The  eddy  dif- 
fusivity  is  similarly  defined  such  that 


-a  c .  , 
P  ij 


3T 

8x 


where  is  the  turbulent  Prandtl  number.  It  should  be  noted  that 

Spraggs  did  not  allow  for  the  effect  of  stratification,  through  the 
Richardson  number,  on  the  eddy  coefficients  in  his  initial  work,  but 
did  indicate  that  such  a  modification  would  be  made  later. 

138.  The  computational  grid  employed  is  one  such  that  the 
velocity  components  are  defined  on  the  cell  faces;  whereas,  the 
thermodynamic  variables  are  defined  at  the  cell  center.  Thus,  the 
grid  is  in  essence  a  grid  similar  to  that  employed  by  Leendertse  (1967). 

139 •  Boundary  conditions  at  solid  walls  are  treated  as  no-slip. 
Thus,  the  normal  velocity  at  a  wall  is  set  to  zero,  and  its  value  at 
one  grid  point  outside  the  wall  is  set  as  the  negative  of  its  value 
at  the  first  interior  point.  Tangential  velocities  are  not  defined 
at  the  wall.  However,  in  order  to  model  the  effect  of  a  no-slip 
wall,  its  value  at  one  grid  point  outside  the  wall  is  taken  to  be 
the  negative  of  its  value  at  one  grid  point  inside.  Both  inflow 
and  outflow  boundaries  are  assumed  to  be  forced.  At  the  surface, 
velocities  are  set  using  a  wind  stress  condition.  The  temperature 
field  at  all  boundaries  except  the  free  surface  is  assumed  to  have 
a  zero  gradient;  whereas,  surface  temperatures,  of  course,  are 
determined  from  the  surface  hent  exchange  determined  by  prevailing 
atmospheric  conditions. 

l40.  In  the  solution  of  the  velocity  and  temperature  fields, 
forward  differences  are  used  to  replace  time  derivatives.  Roache's 
second  upwind  differencing  scheme  is  used  to  replace  the  advective 
terms.  Thus  (see  Figure  5), 


Vr  :  Vl 

Ax 


where  aD 

r\ 


and 


depend  on  the  sign  of  the  convecting  velocities. 
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Centered  differences  are  employed  in  the  representation  of  the  diffu¬ 


sive  terms.  Therefore,  the  finite  difference  scheme  is  in  essence 
2 

0(At,Ax  )  .  Various  schemes  of  higher  order  were  investigated  by 
Spraggs  and  Street,  e.g.,  the  leapfrog,  the  Adams-Bashford,  and  Fromm's 
(see  Roache  1972)  second  order  schemes.  The  leapfrog  scheme  was  dis¬ 
carded  because  of  the  time-splitting  nature  of  the  solution,  while 
Fromm's  method  was  not  used  due  to  the  large  percentage  of  boundary 
cells  encountered  in  3-B  modeling  where  the  method  uses  centered  spatial 
differencing.  Such  a  scheme  was  found  to  be  unacceptable  near  bound¬ 
aries  with  large  forced  outflows.  A  similar  conclusion  was  arrived  at 
during  computer  experimentation  with  the  2-D  Edinger  and  Buchak  (1979) 
model  (page  %).  Spraggs  and  Street  indicate  that  the  necessary  coding 
for  the  Adams-Bashford  method  remains  in  the  basic  numerical  model  for 
future  development  and  testing. 

Eraslan's  3-D  Discrete  Element  Model 


lUl .  Eraslan*  is  currently  working  on  a  fully  three-dimensional 
heat-conducting  hydrodynamic  model  for  the  Oak  Ridge  National  Labora¬ 
tory.  The  code  will  be  a  fully  convective  model  with  the  complete 
vertical  momentum  equation  retained.  The  basic  solution  technique  will 
employ  an  explicit  time-differencing  scheme  along  with  the  previously 
discussed  concept  of  discrete  elements.  Therefore,  his  formulation  will 
employ  integral  forms  of  the  governing  conservation  equations  applied  to 
variable-sized  discrete  elements  that  span  user-specified  flow  regions. 
At  the  present  time,  there  is  no  published  information  on  the  develop¬ 
ment  of  the  model. 


Blumberg  and  Mellor's  3-D  Model 


lU2.  After  the  initial  writing  of  this  report,  a  three- 
dimensional  heat-conducting  coastal  model  developed  by  Blumberg  and 


*  Personal  communication.  May  1979,  Arsev  Eraslan,  Chief  Scientist, 
Hennington,  Durham,  and  Richardson,  Knoxville,  Tenn. 
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Mellor  (1979)  at  rrinceton  University  was  brought  to  the  at t  :  r. 

the  author.  It  appears  that  the  model  is  in  an  early  .  tare  ,-f  a;  ;  .  ;  ■  - 
tion,  with  only  preliminary  tests  in  the  Gulf  01'  Mexico  savin  •. 
made . 

IU3.  The  basic  equations  solved  are  statements  >f  the  •t:..a-rv;- 
tion  of  fluid  mass,  momentum,  and  energy  along  with  the  conservat.  f 

salt  equation.  The  energy  equation  is  written  in  term.;  of  tempera4 are , 
and  thus  the  equation  of  state  relates  the  fluid  density  to  lot:,  tem¬ 
perature  and  salinity.  The  basic  Boussinesq  and  hydrostatic  pressure 
assumptions  are  made. 

144.  The  model  employs  two  concepts  previously  discussed  in  con¬ 
nection  with  the  Simons  and  Lick  models.  Similar  to  the  Simons  model, 
the  external  flow  is  computed  separately  from  the  internal  flow.  The 
external  mode,  an  essentially  two-dimensional  calculation,  requires  a 
short  integrating  time  step;  whereas,  the  three-dimensional,  internal 
mode  can  be  executed  with  a  long  step.  The  result  is  a  fully  three- 
dimensional  code  that  includes  a  free  surface.  Similar  to  the  Lick 
model,  the  vertical  coordinate  is  transformed  into  a  a  coordinate 
system  with  20  levels  in  the  vertical.  The  model  developers  state, 

"With  such  a  transformation,  the  environmentally  important  continental 
shelf,  shelf  break,  and  slope  can  be  well  resolved."  Furthermore,  the 
model  allows  for  variable  grid  spacing  in  the  o  coordinate  for  in¬ 
creased  resolution  in  the  surface  and  bottom  layers. 

ll*5.  Rather  than  employing  the  same  concept  of  eddy  coefficients 
as  utilized  by  all  the  other  models  investigated,  a  second  moment  model 
of  small-scale  turbulence  as  developed  by  Mellor  and  Yamada  (1977)  is 
employed.  Diffusive-type  terms  proportional  to  second  derivatives  in 
the  basic  equations  are  retained  only  in  the  vertical  direction.  The 
developers  indicate  that  they  believe  relatively  fine  vertical  resolu¬ 
tion  results  in  a  reduced  need  for  horizontal  diffusion;  i.e.,  horizon¬ 
tal  advection  followed  by  vertical  mixing  effectively  acts  as  a  horizon¬ 
tal  diffusion  in  a  real  physical  sense. 

1^6.  At  the  surface,  the  wind  stress,  net  heat  flux,  and  net 
evaporation-precipitation  freshwater  flux  are  accounted  for.  Bottom 


boundary  conditions  on  the  velocity  components  are  supplied  by  matcain.-- 
the  solution  to  the  logaritlimic  law  of  the  wall. 

ll+7.  Time  differencing  is  the  conventional  leapfrog  technique. 
However,  the  scheme  is  quasi-implicit ,  since  the  vertical  diffusive 
terms  are  evaluated  at  the  forward  time  level.  Thus,  small  vertical 
spacing  is  permissible  near  the  surface  without  the  need  to  reduce  the 
time  increment  or  restrict  the  magnitude  of  the  mixing  coefficients. 

The  spatial  differencing  is  not  discussed,  but  Blumberg  and  Mellor 
1,1979)  state  that  the  overall  solution  is  accurate  to  the  second  order 
in  space  and  time. 

lU8.  As  previously  discussed,  leapfrog  time  differencing  intro¬ 
duces  a  tendency  for  the  solutions  at  even  and  odd  time  liner,  to  split. 
The  time-splitting  here  is  removed  by  the  use  of  a  weak  filter  where 
the  solution  is  smoothed  at  each  time  step  by 


F 


n 

s 


2Fn  + 


where  a  =  1/10  and  F  is  a  smoothed  solution.  This  technique  intro- 

s 

duces  less  damping  than  either  the  Euler  backward  or  forward  stepping 
techniques  (see  Roache  1972). 
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.'’ART  V:  TWO-DIMENSIONAL  VERTICAL  FLOW  HYDRODYNAMIC  MODELS 


IL9.  Various  two-dimensional  numerical  hydrodynamic  models  have 
been  studied  and  range  from  vertically  averaged  models  to  laterally 
averaged  models  to  pure  two-dimensional  vertical- longitudinal  models  in 
which  the  width  is  constant.  There  are  many  two-dimensional  vertically 
averaged  estuarine  models  in  existence.  However,  since  such  models  are 
not  applicable  to  stratified  reservoir  flows  in  which  the  allowance  for 
a  variable  density  and  variations  in  the  vertical  direction  are  crucial, 
they  are  not  discussed  here.  The  only  interest  in  such  models  in  conneo 
tion  with  this  study  was  in  the  numerical  techniques  employed.  All  of 
the  models  studied,  except  for  one  that  employs  the  finite  element 
method ,  utilize  the  finite  difference  method  for  solving  the  governing 
f-D  equations.  Vnlike  all  of  the  3-D  models,  which  were  explicit  models 
some  of  the  1-1  models  employ  an  implicit  or  perhaps  semi- implicit  time 
integration  scheme  so  that  time  steps  much  larger  than  that  -iven  by 
the  dour ant  condition  are  allowed.  A  few  of  the  2-Z  vertical  models  in¬ 
vestigated  were  developed  originally  for  application  to  reservoirs,  and 
thus  surface  heat  exchange  and  the  variability  of  density  with  tempera¬ 
ture  are  treated.  The  Edinger  and  Buchak  (1979),  ’Waldrop  and  Farmer 
(197b),  and  Roberts  and  Street  (1975)  models  are  examples.  fther  .7-: 
vertical  models,  such  as  those  of  Hamilton  (1975)  and  clumbers  (  1.  '<7* ) 
were  initially  developed  for  salinity-stratified  estuaries  and  ndli- 
tional  modifications  would  bo  needed  for  application  to  reservoirs.  Ad¬ 
ditional  density-varying  models  that  consider  flow  in  a  vortical  : lane 
have  been  investigate!  and  include  those  of  Thompson* ,  • csoi  ion,  Inc.,** 
Norton,  King,  and  Oriob  (i  >73),  and  .'lot ta  ot  nl.  ( 1 9i” .* ' . 

Hamilton's  Estuary  Model 

150.  Hamilton's.  (1975)  -0  model  war  love  lore  i  to  represent  t.i.e 

*  Personal  communication,  April  lpi'o,  .  y.  Phompson,  Missis.;;;  i 
State  University,  Mississippi  State,  Miss. 

**  Personal  communication.  May  197-?,  •'•••rsonnel  ,:>f  :  'ail  f 


vert  icul  .•ir’u't of  current,  and  salinity  along  an  estuary  of  varying 
wi  ;ih  and  depth  but  with  a  reet angular  cross  section,  i.e.,  B  4  • 

A  feature  of  the  model  is  the  numerical  approach  to  the  basic  equations, 
which  considers  the  depth-dependent  variables  as  continuous.  This  is 
in  contrast  with  layered  models  where  the  equations  are  integrate :  v-  r 
tile  se:  urate  layers  and  exc'aans.e  of  momentum  and  salt  between  lay.-;  s 
i  urometorizcd  In  terms  of  the  ::.ean  velocities  and  salinities  of  Ins¬ 
incere.  dan.i  Lton  indicates  tiiat  a  continuum  approach  allows  Letter 
treat::. ei. t  of  the  surface  and  botto::.  boundary  conditions. 

iti.  ihe  i-aaic  3-1  equations  are  reduced  to  a  set  of  laterally 
averaged  il-i.  equations  as  previously  outlined.  Tlic  only  difference  here 
is  taat  the  width  is  not  a  function  of  the  vertical  coordinate  and, 
thus,  derivatives  of  the  width  with  respect  to  the  ver' ..cal  coordinate 
z  are  zero,  Tne  resulting  laterally  averaged  equations,  with  the 
sou ins inesq  approximation  and  the  hydrostatic  pressure  assumption,  are 
written  in  nonconservative  form.  Salinity  is  related  to  the  density 
through  a  linear  equation  of  state. 

152.  Boundary  conditions  at  the  head  of  freshwater  flow  consist 
of  a  vertical  velocity  profile  and  zero  salinity.  At  the  ocean  boundary 
the  tidal  elevation  is  prescribed  as  a  function  of  time,  and  the  salinit 
is  specified  to  be  that  of  trie  ocean.  If  does  not  appear  that  Hamilton 
delineates  an  inflow  and  an  outflow  boundary  at  the  ocean  end.  To  con¬ 
serve  salt,  the  vertical  salinity  gradient  at  the  estuary  surface  and 
the  bed  is  set  to  zero.  Surf  ice  wind  stress  is  neglected,  and  the  bet¬ 
ter.  stress  is  assumed  to  obey  the  quadratic  friction  law  such  that  the 

; . atom  stress  is  related  to  the  velocity  at  a  distance  above  the  bottom 
representative  of  the  frictional  layer,  e. ~. ,  1  m. 

153.  as  in  all  hydrostatic  models,  the  vertical  component  of 
velocity  is  obtained  Vy  solving,  the  laterally  averaged  incompressibility 
condition  from  the  bottom  upward.  The  equation  for  the  free  surface  is 
obtained  by  vertically  integrating  the  incompressibility  equation.  One 
restriction  imposed  on  the  free  surface  by  the  code  logic  is  that  the 
snrfa:e  elevation  does  not  differ  by  more  than  the  vertical  grid  spacing 
(assumed  to  be  constant)  between  successive  horizontal  grid  points. 


15*+.  The  finite  difference  grid  is  such  that  salinity,  the  verti¬ 
cal  velocity,  and  the  vertical  eddy  viscosity  and  diffusivity  are  de¬ 
fined  at  the  center  of  a  cell;  whereas,  the  horizontal  velocity  is  de¬ 
fined  at  the  cell  corners.  This  is  illustrated  below. 


155.  The  time  integration  is  a  combination  forward  and  backward 
time-differencing  scheme  such  that  the  diffusive  and  frictional  terms 
in  the  conservation  of  salt  and  momentum  equations,  respectively,  are 
taken  at  the  n  +  1  time  level,  while  all  other  terms  such  as  the 
advective  terms  are  taken  at  the  n  time  level.  Spatial  differences 
are  replaced  by  centered  differences,  except  in  the  horizontal  advective 
term  of  the  conservation  of  salt  equation,  i.e.,  uSs/'+x  ,  in  which 
Hamilton  appears  to  make  use  of  fioache's  (1972)  first,  upwind  differencing. 
Thus,  the  finite  difference  scheme  is  in  essence  0(At ,Ax)  .  The  basic 
stability  criterion  is  the  Courant  condition.  Therefore ,  even  though 
the  scheme  might  be  called  a  semi-implicit  one  because  the  second  deriva¬ 
tive  terms  are  taken  at  the  n  +  1  time  level,  which  does  remove 
diffusive-type  stability  criteria,  the  scheme  proh.nl ly  offers  ;.e  real 
stability  advantages  over  a  purely  explicit  scheme. 


dij 


Blumberg's  2-D  Laterally  Averaged  Estuary  Model 

156.  Like  the  Hamilton  model,  Blumberg's  (1975)  laterally  averaged 
model  was  developed  for  application  to  an  estuary.  Thus,  the  density 
is  related  to  the  salinity  through  an  equation  of  state  and,  of  course, 
no  surface  heat  exchange  is  included,  since  temperature  is  not  modeled. 
Unlike  the  Hamilton  model,  however,  this  model  does  not  assume  a 
rectangular  cross  section  and  thus  B  =  B(x,z)  . 

157-  Additional  assumptions  made  to  the  basic  equations,  which 
are  written  in  conservative  form,  are  that  the  pressure  is  hydrostatic, 
the  Boussinesq  approximation  is  applicable,  and  that  eddy  coefficients 
can  be  employed  to  represent  the  effect  of  subgrid-scale  motions.  Verti¬ 
cal  velocities  are  thus  computed  from  the  incompressibility  condition, 
and  the  free  surface  equation  re '"Its  from  a  vertical  integration  of 
the  equation  for  incompressibility. 

158.  Boundary  conditions  imposed  consist  of  the  inflow  of  fresh 
water  with  zero  salinity  at  the  head  of  the  estuary,  whereas,  salinity 
and  tidal  elevations  are  specified  at  the  ocean  end.  Unlike  Hamilton, 
Blumberg  allows  for  the  ocean  boundary  to  be  alternately  an  inflow  and 
then  an  outflow  boundary.  When  inflow  occurs,  the  salinity  is  set  to 
be  that  of  the  ocean;  during  outflow,  it  is  determined  from  an  extrapola¬ 
tion  of  values  inside.  To  prohibit  the  flux  of  salt  through  the  surface 
and  the  bottom,  the  vertical  salinity  gradients  are  set  to  zero  at  those 
locations.  The  boundary  condition  on  the  velocity  at  the  surface  is 
determined  from  the  wind  stress.  Similarly,  the  bottom  stress  determines 
the  boundary  condition  at  the  bottom.  Extrapolation  from  the  hydraulic 
theory  of  flow  in  open  channels  allows  the  friction  acting  on  a  tidal 
current,  because  of  the  estuary's  bottom,  to  be  expressed  using  the 
quadratic  law 

T  =  ku | u | 

where  u  is  evaluated  1  m  away  and  k  depends  primarily  on  the 
boundary  roughness. 

-T 


159.  The  basic  finite  difference  grid  is  of  the  MAC-type  need  by 
Leendertse  (1967).  Pressure  and  salinity  are  defined  at  cell  centers, 
but  velocities  are  defined  on  the  faces  of  the  cells.  Water  surface 
elevations  are  defined  on  columns  corresponding  to  cell  centers.  Dira- 
ilar  to  the  layered  approach  of  Leendertse,  trie  governing  equations  are 
integrated  vertically  over  each  layer-  where  the  thickness  of  each  layer 
is  constant  except  for  the  top  one.  The  top  layer,  of  course,  contains 
the  influence  of  the  surface  gravity  wave,  and  its  thickness  varies  in 
time  and  space. 

160.  Both  a  horizontal  and  a  vertical  eddy  viscosity  coefficient 
as  well  as  a  horizontal  and  vertical  diffusivity  coefficient  are  com¬ 
puted.  The  horizontal  coefficients  are  computed  from 


where 

C  =  adjustable  constant 

while  the  vertical  coefficients  are  related  to  the  Richardson  number  in 
the  following  manner: 


and 

A 

v 

EV  "  1  +  R. 

1 

where  is  the  eddy  diffusivity  and  is  the  eddy  viscosity,  k 

is  a  constant  whose  value  is  ~0.10  and  is  a  critical  Richardson 

c 

number  taken  to  be  10.  It  should  be  remembered  that  Blumberg's  model 
was  developed  for  an  estuary.  Therefore,  the  functional  form  of  the 
coefficients  above  are  probably  not  applicable  to  deep  reservoirs.  As 
was  done  in  the  3-D  quasi-static  models,  the  eddy  diffusivity  is  assumed 
large  when  unstable  stratification  develops.  The  salinity  in  the  un¬ 
stable  layers  is  replaced  by  the  averaged  value  of  the  adjacent  layers. 

l6l.  The  time-integration  scheme  is  a  centered  difference  or 
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leapfrog  scheme,  except  for  the  diffusive  and  frictional  terms,  which 
are  taken  at  the  old  time  step.  All  spatial  derivatives  are  replaced  by- 
centered  differences.  Thus,  as  with  the  Simons  3-D  model,  the  finite 
difference  scheme  is  almost  0(At  ,Ax  )  .  As  previously  noted,  the  use 
of  centered  differences  in  time  and  space  results  in  a  second  order  dif¬ 
ference  equation  as  the  approximation  to  a  first  order  differential 
equation,  and  the  solutions  at  odd  and  even  time  lines  tend  to  split. 
Blumberg  attempts  to  remove  this  time-splitting  through  averaging  re¬ 
sults  from  three  successive  time  steps  with  weights  of  0.25,  0.50,  and 
0.25,  respectively,  every  25  time  steps. 

162.  The  centered  difference  time-integration  scheme  has  the 
property  of  not  introducing  artificial  horizontal  diffusion  and  vis¬ 
cosity.  Thus,  to  control  nonlinear  instabilities,  damping  must  be  input 
into  the  scheme.  This  is  the  major  reason  for  incorporating  the  ex¬ 
pressions  previously  given  for  the  horizontal  diffusivity  and  viscosity, 
A^  and  ,  respectively. 

Poseidon's  2-D  Vorticity-Stream  Function  Model 

163.  There  are  no  published  reports  on  Poseidon's*  2-D,  longi¬ 
tudinally  and  vertically  dimensional,  variable  density  model.  The  major 
reasons  for  noting  the  model's  existence  are  first  because  it  is  the 
only  hydrodynamic  model  discovered  that  is  based  on  the  vorticity-stream 
function  representation  of  the  governing  equations  and  secondly,  because 
of  the  manner  in  which  the  advection  terras, 

3(ug)  +  3(vg) 

3x  3y 

where  <;  is  vorticity,  are  numerically  modeled.  As  noted  before,  the 
oasic  problem  with  these  terms  is  that  of  achieving  numerical  stability 
without  numerical  diffusion.  The  Poseidon  code  uses  a  flux-corrected 
transport  algorithm  called  SHASTA  (Sharp  and  Smooth  Transport  Algorithm). 

*  Personal  communication,  May  1978,  Personnel  of  Poseidon,  Inc.,  Calif. 
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Fluxes  are  first  advected  according  to  a  scheme  that  is  stable  but  dif¬ 
fusive,  e.g.,  the  two-step  Lax-Wendroff  algorithm.  Then  the  amount  of 
numerical  diffusion  is  computed  at  each  grid  point,  and  the  appropriate 
amount  of  antidiffusion  flux  is  applied  to  each  cell,  provided  no  new 
extrema  are  created.  A  discussion  of  SHASTA  is  given  by  Boris  and 
Book  (1973).  Again  it  should  be  noted  that  such  a  model  would  not  be 
applicable  to  a  reservoir  containing  multiple  outlets . 

Slotta  et  al.’s  2-D  NUMAC  Model 


l6U.  A  group  directed  by  Slotta  ( Slotta  et  al.  1969)  at  Oregon 
State  University  has  developed  the  computer  model  NUMAC  (Nonhomogeneous 
Unconfined  Marker  and  Cell)  for  analyzing  transient,  incompressible, 
variable  density,  viscous  flows  with  a  free  surface.  As  the  name  im¬ 
plies,  the  model  is  based  upon  the  MAC  method  developed  by  Welch  et  al. 
(1966),  which  uses  a  mixed  Euler ian-Lagrangian  scheme.  In  this  scheme, 
the  velocity  and  pressure  are  considered  as  Eulerian  variables  defined 
at  the  mesh  points  of  a  fixed  grid,  but  the  density  is  considered  a 
Lagrangian  variable  localized  to  fluid  particles.  It  appears  that  the 
major  differences  between  NUMAC  and  MAC  lie  in  NUMAC' s  ability  to  better 
handle  inlets  and  outlets  and  in  the  use  of  the  SOR  technique  for 
solving  the  Poisson  equation  for  the  pressure. 

165.  The  basic  Navier-Stokes  equations  for  laminar  flow  written 
in  the  vertical  and  longitudinal  directions,  in  the  conservative  form, 
are  solved  along  with  the  conservation  of  mass  equation.  The  Boussinesq 
approximation  is  not  made,  and  thus,  the  density  is  actually  solved  for 
from  a  transport  equation  with  p  as  the  dependent  variable.  However, 
the  incompressibility  condition  is  still  invoked  in  the  derivation  of 
the  Poisson  equation  for  the  pressure. 

166.  Many  different  types  of  boundary  conditions  are  allowed. 

At  material  boundaries,  the  normal  component  of  the  velocity  vanishes . 

At  a  free  surface,  the  boundary  conditions  are  that  the  normal  and 
tangential  components  of  the  stress  must  vanish.  Two  inlet  velocity 
boundary  conditions  are  allowed.  One  holds  the  inlet  velocity  constant. 


while  the  other  requires  the  normal,  derivative  to  vanish.  The  normal 
derivative  of  the  density  at  an  outlet  is  set  to  zero.  Both  slip  and 
no-slip  solid  boundaries  are  allowed  with  the  derivative  boundary  con¬ 
dition  on  the  pressure  determined  from  the  momentum  equation. 

167.  The  finite  difference  scheme  is  basically  one  in  which  the 
time  derivatives  are  replaced  by  forward  differences  and  the  spatial 
derivatives  by  centered  differences.  However,  it  does  appear  that 
Roache's  (1972)  second  windward- type  differencing  is  used  in  the  eval¬ 
uation  of  momentum  flux  terms  such  as  3(puv)/3y  ,  etc.  Thus,  theoret- 
ically,  the  scheme  is  close  to  0(At,  Ax  )  . 

168.  As  noted  previously,  the  MAC  calculations  are  a  combination 
of  Eulerian  and  Lagrangian  steps.  The  NUMAC  computation  cycle  is  sum¬ 
marized  in  the  following  steps: 

a.  Compute  new  densities  from  the  mass  transport  equation. 

b_.  Using  new  densities,  solve  Poisson  equation  for  the 
pressure . 

£.  Using  new  densities  and  pressures,  calculate  new  veloc¬ 
ities  from  momentum  equations . 

d.  Move  the  Lagrangian  particles  by  use  of  the  new 
velocities . 

e_.  Calculate  new  densities  and  viscosities  at  the  mesh 
points  by  averaging  the  densities  and  viscosities  of 
the  particles  that  now  surround  each  mesh  point. 

f_.  Compare  this  density  with  the  value  computed  in  step  a. 
If  different,  go  to  step  b_  with  these  densities.  If 
they  are  essentially  the  same,  continue. 

£.  Recompute  the  pressure  from  the  Poisson  equation. 

h_.  Recompute  the  velocities  from  the  momentum  equations. 

£.  Move  the  particles  using  velocities  from  step  h_. 

j_.  Increment  the  time  and  go  to  step  a. 

169.  Several  stability  criteria  for  this  explicit,  scheme  are  pre¬ 
sented;  however,  once  again  the  basic  criterion  is  related  to  the  speed 
of  a  gravity  wave. 

170.  Obviously,  NUMAC,  or  any  of  the  related  MAC  codes,  is  an 
extremely  powerful  numerical  model  for  analyzing  variable  density  fluid 
flows,  since  the  model  is  fully  convective.  However,  computing  times 

89 


required  for  long-term  transient  problems  are  excessive  due  to  the 
explicit  time  differencing  plus  solving  a  Poisson  equation  for  the 
pressure.  Slotta  indicates  that  one  time  cycle  on  a  (Control  Data 
Corporation)  CDC  6600  computer  requires  7  sec  for  a  problem  with 
800  cells  and  3000  particles. 

Norton,  King,  and  Orlob's  2-D  Vertical 
Flow  FEM  Model— RMA-7 

171.  Under  a  contract  with  the  Walla  Walla  District  of  the  U.  S. 
Army  Corps  of  Engineers,  Water  Resources  Engineers,  with  Norton  and  King 
as  principal  investigators,  developed  two  2-D  hydrodynamic  models  using 
the  finite  element  method  for  obtaining  numerical  solutions  of  the  gov¬ 
erning  flow  equations  (Norton,  King,  and  Orlob  1973).  One  of  the  models 
is  a  variable  density,  laterally  averaged  model  that  describes  the  be¬ 
havior  of  velocity,  temperature,  and  pressure  in  the  vertical  plane. 

172.  The  basic  equations  solved  are  the  2-D  laterally  averaged 
horizontal  and  vertical  momentum  equations  along  with  the  continuity 
equation  reduced  to  the  incompressibility  condition  as  a  result  of  the 
Boussinesq  approximation  and  an  energy  equation  written  in  terms  of 
temperatures.  These  four  equations  along  with  an  equation  of  state 
relating  the  fluid  density  to  the  temperature  are  solved  for  the  five 
unknowns — u  ,  v  ,  T  ,  P  ,  and  p  . 

173.  The  exchange  of  energy  to  the  unresolvable  turbulent  eddies 
is  accomplished  through  the  use  of  eddy  coefficients,  which  are  treated 
as  constants  within  each  element  but  can  vary  from  element  to  element. 

It  should  be  noted  that  unlike  most  models,  the  off-diagonal  terms  of 
the  eddy  viscosity  tensor  are  retained. 

17^.  The  equations  are  written  in  the  nonconservative  form  with¬ 
out  the  usual  hydrostatic  approximation.  Thus ,  the  complete  vertical 
momentum  equation  is  retained  and  the  model  is  a  fully  convective  model. 

175.  The  governing  equations  are  solved  by  the  finite  element 
method  using  Galerkin's  method  of  weighted  residuals.  A  mixed  set  of 
basic  functions  is  employed  in  the  overall  permutation.  ..uadratio 
function.:  are  used  for  all  state  variables  except  pressure  where  a 
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linear  function  is  used.  The  linear  pressure  function  implies  a  con¬ 
stant  element  density,  which  is  calculated  as  a  function  of  average 
nodal  temperatures.  An  implicit,  Newton- Raphson  computation  scheme  is 
employed  to  achieve  a  solution  to  the  set  of  nonlinear  equations  that 
define  the  model.  The  resulting  computer  program  accommodates  tri¬ 
angular  and/or  quadrilateral  isoparametric  elements. 

176.  Both  a  bottom  stress  term  and  a  wind  shear  term  are  incor¬ 
porated  in  tile  bottom  and  top  row  of  elements,  respectively.  The  use 
of  the  isoparametric  formulation  with  interelement  geometric  slope  con¬ 
tinuity  allows  the  user  to  specify  slip  or  parallel  boundary  flows.  In 
addition,  no-slip  walls  can  be  easily  handled  since  zero  values  of  u 
and  v  would  be  inserted  at  the  proper  nodes  of  boundary  elements.  The 
surface  heat  flux  at  the  air-water  interface  is  computed  through  the  use 
of  tiie  coefficient  of  surface  heat  exchange  and  local  equilibrium  tem¬ 
perature  as  calculated  from  meteorological  data. 

177.  A  recent  version  of  the  model  accounts  for  the  movement  of 

the  free  surface,  alt:  ough  in  a  very  limited  fashion,  since  the  movement 
must  be  stipulated  by  the  user.  The  free  surface  pressure  boundary  con¬ 
dition  is  based  upon  the  assumption  of  a  locally  surface  so  that 

the  pressure  boundary  condition  is  for  atmospheric  pressure.  The  model 
developers  are  currently  incorporating  into  the  model  a  procedure  for 
internally  computing  the  location  of  the  free  surface  utilizing  the 
atmospheric  pressure  boundary  condition. 

178.  As  noted  in  previous  discussions,  finite  element  models  for 
transient  problems  require  large  computing  times.  Therefore,  such 
models  may  not  be  applicable  to  the  simulation  of  the  natural  stratifi¬ 
cation  cycle  of  a  reservoir  for  economic  reasons. 

179-  It  might  be  noted  that  although  laterally  averaged  models 

provide  a  better  representation  of  real  reservoirs  than  pure  2-B  models, 

the  momentum  flux  through  an  outlet  at  the  dam  is  not  accurately  modeled. 

In  the  horizontal  momentum  equation,  the  horizontal  advection  of  mortien- 

-2  - 

turn  is  represented  by  pq3(u  B)/3x  ,  where  u  is  the  laterally  averaged 

velocity  in  the  x  direction.  In  actuality,  the  momentum  flux  passing 
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a  cross  section  is  not  pQU  B  ,  but  instead  is  given  by  the  integral 


a  value  of  1.0  is  reached.  However,  the  usual  procedure  is  to  assume  a 
value  of  1.0  everywhere,  in  which  case  the  flux  of  momentum  at  the  down¬ 
stream  boundary  is  not  properly  modeled. 

Thompson's  2-D  Model — VESSEL 

180.  This  is  a  laterally  averaged  2-D  model  that  is  currently 
being  developed  to  assist  the  Corps  in  selective  withdrawal  studies. 
Because  of  the  concern  for  the  quality  of  water  downstream  of  reser¬ 
voirs,  there  is  a  growing  effort  to  control  the  quality  of  water  re¬ 
leased  from  reservoirs.  The  concept  of  controlling  the  quality  released 
from  a  density-stratified  impoundment  is  called  "selective  withdrawal." 
Because  the  quality  of  water  and  its  density  can  vary  from  the  surface 
to  the  bottom  of  a  lake,  it  is  often  possible  to  selectively  withdraw 
the  most  desirable  qualities.  A  basic  problem  is  to  determine  before 
construction  whether  the  design  of  an  outlet  will  provide  the  desired 
selective  withdrawal  characteristics. 

181.  An  empirical  method  developed  by  Bohan  and  Grace  ( 1969 )  can 
be  utilized  for  selective  withdrawal  predictions  for  simplified  outlet 
and  approach  geometries.  However,  for  complex  geometries,  physical 
and/or  mathematical  models  are  required. 

182.  Thompson's  model  utilizes  the  concept  of  boundary-fitted 
coordinates  to  obtain  a  solution  of  the  governing  flow  equations  on  a 
nonorthogonal  curvilinear  coordinate  system.  The  coordinate  system  is 
generated  from  the  elliptic  generating  system 


S  +  £  =  P 

xx  yy 


n  +  h  =  Q 
xx  yy 


where  P  and  Q  are  functions  chosen  to  cause  the  ^  ,  q  coordinate 


lines  to  concentrate  as  desired.  With  one  coordinate  being  specified 
as  constant  on  the  boundaries  and  a  distribution  of  the  other  specified, 
a  coordinate  system  that  follows  all  boundaries,  no  matter  how 


irregular,  results.  A  rather  detailed  disc  us.,  ion  of  method  an  1  it.' 
v oss i bio  application  to  hydro  lynamic  problem.:  is  presence  i  :  y  Pc  .nr:  •,  :.n 
arid  Thompson  ( 1973 ) .  in  addition,  an  extent  ivo  list  of  references  in¬ 
scribing  Thompson '  s  work  with  the  technique  is  presentee  is.  the  work 
cited . 

193.  Tiie  next  step  is  the  development  of  a  numerical  model  to 
solve  the  governing^  fluid  flow  equations  on  the  coordinate  system  com¬ 
puted  above.  Such  a  .model  will  be  able  to  accurately  model  the  in¬ 
fluence  of  boundary  geometry  on  the  developing  flow. 

l8U.  The  basic  laterally  averaged  2-D  equations  solved  in 
Thompson's  model  are  the  iiavier-Stokes  equations,  mass  conservation, 
energy  conservation,  and  an  equation  of  state  relating  temperature  and 
density.  These  equations  are  transformed  to  the  £  ,  n  system  in  a 
fully  geometrically  conservative  form  sucii  that  the  finite  difference 
representation  is  equivalent  to  the  discrete  element  method.  Essen¬ 
tially  nc  assumptions  other  than  assuming  an  incompressible  fluid  are 
applied  to  the  basic  equations;  e.g.  ,  the  Boussinesq  approximation  is 
not  made  and  the  model  is  fully  convective  with  the  vertical  velocity 
obtained  from  the  full  vertical  momentum  equation.  In  the  vicinity  of 
outlets,  vertical  accelerations  may  become  large  and  a  solution  of  the 
full  vertical  momentum  equation  is  probably  required. 

185.  The  pressure  is  computed  using  Chorin's  method.  This  method 
is  based  upon  the  concept  that  if  a  fluid  is  incompressible,  the  func¬ 
tion  of  the  pressure  is  to  insure  that  the  velocity  field  satisfies  the 
incompressibility  condition,  i.e.,  V  •  v  =  0  .  An  iterative  algorithm 
for  the  pressure  field  is  thus  set  up  such  that 
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where  the  pressure  field  at  time  step  n  in  determined  sucii  that  the 
velocity  field  at  time  step  n  +  1  will  satisfy  incompressibility.  The 
advantage  of  Chorin's  method  over  the  use  of  a  Poisson  equation  is  that 
only  velocities  are  required  on  the  boundaries,  rather  than  pressure 
and/or  velocity  derivatives. 
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186.  The  finite  difference  grid  is  such  that  both  velocity  com- 
ponents  are  defined  at  cell  corners,  while  the  pressure,  temperature ,  and 
density  are  defined  at  the  center  of  a  cell.  It  should  be  noted  that 
the  transformation  of  the  equations  into  the  boundary- fitted  coordinate 
system  is  such  that  all  computations  are  performed  on  a  rectangular 
(C  ,  h)  grid  with  square  grid  mesh. 

137.  The  model  being  developed  will  be  extremely  general  so  that 
any  number  of  inlets  and/or  outlets  can  lie  on  any  boundary.  In  addi¬ 
tion,  any  number  of  bodies  can  lie  in  the  interior  of  the  field,  with  a 
constant  coordinate  line  following  each  body.  Boundary  conditions  can 
be  either  slip  or  no-slip  on  solid  boundaries,  with  the  option  of  either 
specifying  wall  temperatures  or  the  heat  transfer  rate  at  such 
boundaries . 

188.  The  bas  ic  finite  difference  scheme  utilizes  second  order 
backward  differences  to  replace  time  derivatives  and  centered  differ¬ 
ences  to  replace  spatial  derivatives.  The  model  will  allow  the  option, 
however,  of  selecting  windward  differencing  of  advective  terms.  The 

finite  difference  scheme  is  thus  fully  implicit  and  of  0(At  ,  AC  )  or 
2  2 

almost  0(At  ,  AC  )  ,  depending  upon  whether  Roache's  (1972)  first  or 
second  differencing  is  employed.  The  SOR  iterative  method  with  a  vari¬ 
able  optimum  acceleration  parameter  field  is  utilized  to  obtain  a 
solution. 

189.  With  such  an  unsteady,  fully  convective,  variable  density, 
free  surface  model  that  models  the  flow  phenomena  in  a  natural  coordi¬ 
nate  system  that  fits  the  boundaries  of  the  field,  a  wide  range  of  hy¬ 
draulic  phenomena  can  be  accurately  simulated.  However,  due  to  the 
fully  implicit  nature  of  the  solution  and  the  resulting  iterative  solu¬ 
tion  technique,  the  computing  time  required  for  long-term  simulations 
will  probably  be  large.  For  selective  withdrawal  studies  in  which  only 
the  steady-state  solution  is  sought,  the  computing  cost  should  not  be 

a  major  factor. 

Roberts  and  Street's  2-7  Reservoir  Model 

190.  Roberts  and  Street ' s  (1975)  variable  density  model  is  quite 
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similar  to  T:  ra.ys  and  .'tre.t'.-  model  wit.:,  One  ma,’  ex  - 

t i  on::  being  the  iroppin.-*  of  the  lateral  dimension  ar.  i  the  asca::.p;ti  of 
a  .uyirostat ic  pressure  distribution,  eliminating  t.ne  nee.:  for  a  .'era 
pressure  equation  and  the  attendant  costly  solution  procedure.  The 
basic  finite  difference  gri  i ,  solution  technique,  and  eddy  viscosity 
model  are  all  essentially  the  game  as  employed  in  the  3-1  model,  but 
are  now  reduced  to  two  dimensions.  The  model  is  thus  a  pure  2-1 
vertical-longitudinal  model  in  which  a  varying  width  is  not  allowed. 

191.  With  the  hydrostatic  pressure  assumption,  the  vertical 
velocity  is  solved  from  the  condition  of  incompressibility,  and  a  large 
vertical  diffusivity  is  invoked  to  simulate  convective  overturning, 
which  cannot  be  dealt  with  explicitly.  Unlike  some  of  the  nydrostatic 
models  that  integrate  the  incompressibility  equation  over  the  vertical 
to  yield  an  equation  for  the  free  surface,  Roberts  and  street  determine 
the  free  surface  directly  from  the  kinematic  boundary  condition  at  the 
surface.  As  in  the  3-D  Spraggs  and  liroet  model,  an  implicit  solution 
of  the  surface  equation  is  obtained.  Once  again,  however,  because  cf 
the  lack  of  coupling  between  the  velocity  field  and  she  free  surface  at 
time  level  n  +  1  ,  the  Courant  condition  is  still  the  coni  roll  in¬ 
stability  criterion. 

192.  Limited-slip  solid  boundaries  are  assumed.  The  velocity 
orthogonal  to  the  boundary  i.s  set  tc  zero,  but  the  tangential  velocity 
is  defined  by  the  Chezy-Mannin.-r  formula  for  boundary  shear  stress  such 
that  the  proper  velocity  profile  near  the  boundary  can  be  achieved. 
Forced  flow  boundaries,  of  course,  require  the  specification  of  the 
velocity.  At  solid  boundaries ,  temperature  gradients  are  set  to  zero 
to  model  an  adiabatic  wall. 

193-  At  the  free  surface,  th  velocity  boundary  condition  is 
determined  by  the  wind  stress,  and  temperatures  are  determined  by  a 
surface  heat-exchange  equation. 

191*.  The  basic  finite  difference  scheme  for  the  int‘-rnal  :'j  ov 
utilizes  forward  differencing,  in  .ime  an:  centered  iifferencin.-  Iv. 
space,  except  for  the  aivec 1  i.ve  terms,  where  Roach-  V 
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upwind  differencing  is  utilized.  The  overall  solution  is,  thus,  almost 
of  0(At,  Ax2)  . 


Waldrop  and  Farmer's  TVA  2-D  Reservoir  Model 

195-  Waldrop  and  Farmer's  (1976)  model  is  an  explicit  laterally 
averaged  hydrodynamic  model  for  analyzing  flows  in  stratified  reservoirs 
or  long  river  reaches.  The  model  is  designed  to  accommodate  hourly 
changes  in  boundary  conditions  consisting  of  dam  discharges,  tributary 
inflow  conditions,  steam  plant  intake  and  discharge  conditions,  river 
inflow  rates  and  temperatures,  meteorology  and  wind  shear. 

196.  Very  little  detailed  published  material  on  the  model  exists, 
although  Waldrop  and  Walter  Harper  of  TVA  are  currently  in  the  process 
of  writing  such  a  report.  It  should  be  noted  that  Harper  has  been  respon¬ 
sible  for  most  of  the  coding  and  testing  of  the  model;  thus,  the  model 
should  probably  be  called  the  Waldrop-Harper  model.  From  the  limited 
material  available,  it  appears  that  the  nonconservative  form  of  the 
laterally  averaged  fluid  flow  equations  and  the  temperature  transport 
equation,  in  which  the  Boussinesq  approximation  and  the  hydrostatic 
pressure  assumption  have  been  made,  are  solved.  The  effect  of  turbu¬ 
lence  is  included  through  eddy  coefficients,  which  are  modeled  by  using 
a  mixing  length  theory  as  in  Waldrop  and  Tatom's  3-D  model.  The  retard¬ 
ing  effect  of  stratification  upon  vertical  mixing  is  included  by  damping 
the  vertical  eddy  coefficients  as  a  function  of  the  local  Richardson 
number . 

197-  Free  surface  boundary  conditions  on  the  temperature  and 
velocity  are  provided  by  the  specification  of  the  surface  heat  flux  and 
the  wind  shear,  respectively.  The  surface  heat  flux  qg  is  prescribed 
as  a  quadratic  function  of  the  temperature,  given  as 

q  =  a(t)#T2  +  b(t)#T  +  C(t) 
s  s  s 

where  a  ,  b  ,  and  c  are  coefficients  dependent  upon  meteorological 
conditions,  and  Tg  is  the  surface  temperature;  wind  shear  is  given  by 
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TWIND  =  C*^UWIND  "  V 


V3 


where  C  is  a  coefficient,  UyiND  '”'S  the  w^n<*  velocity,  and  Ug  is  the 
surface  water  velocity. 

198.  The  basic  finite  difference  grid  appears  to  be  a  2-D  version 
of  the  3-D  model,  and  thus,  solid  boundaries  are  treated  as  in  that 
model.  In  other  words,  slip  boundary  conditions  are  assumed  at  vertical 
walls;  whereas,  a  limited-slip  condition  is  applied  at  the  bottom  by 
using  a  logarithmic  profile  to  set  the  velocity  near  the  bottom.  With 
such  a  technique,  the  bottom  never  actually  lies  on  a  grid  point. 

199.  With  the  same  basic  finite  difference  scheme  as  employed  in 
the  3-D  model;  i.e.,  a  form  of  forward  differencing  in  time  and  centered 
differencing  in  space,  except  for  Roache's  (1972)  first  windward  differ¬ 
encing  of  advective  terms,  the  basic  scheme  is  probably  of  0(At,  Ax)  . 

Edinger  and  Buchak' s  Laterally  Averaged 
Reservoir  Model — LARM 

200.  Edinger  and  Buchak' s  (1979)  LARM  (Laterally  Averaged  Reser¬ 
voir  Model)  is  a  numerically  efficient  2-D  laterally  averaged  free  surface, 

variable  density,  heat-conducting  model  developed  for  the  Ohio  River  Divi¬ 
sion,  U.  S.  Army  Corps  of  Engineers,  for  use  in  simulating  flows  in 
stratified  reservoirs.  As  noted  by  Edinger  and  Buchak,  "Such  a  model 
is  needed  in  long,  narrow  reservoirs  that  exhibit  density  flow,  epilim- 
netic  wedges  and  titled  isotherms  and  in  deep  power  plant  discharge 
canals  with  bottom  intrusion  of  cold  water  and  backwater  density  wedges 
from  such  discharges  to  rivers." 

201.  In  the  initial  development  of  the  model,  it  was  anticipated 
that  its  primary  use  would  be  for  long-term  simulations  extending  over 
a  natural  stratification  cycle  of  a  reservoir.  Thus,  it  was  deemed 
necessary  to  develop  a  solution  technique  that  would  allow  for  time 
steps  significantly  larger  than  those  imposed  by  the  free  surface 
gravity  wave.  To  allow  this,  finite  difference  techniques  have  been 
employed  to  solve  the  governing  equations  such  that  the  water  surface 


elevations  are  computed  implicitly,  ^he  velocity  components  in  the 
longitudinal  and  vertical  directions  are  then  computed  explicitly,  and 
finally  the  temperature  field  is  computed  implicitly.  The  density  is 
then,  of  course,  computed  from  an  equation  of  state.  Unlike  the  Roberts 
and  Street  (1975)  model,  which  also  implicitly  computes  the  water  sur¬ 
face,  Edinger  and  Buchak's  model  couples  the  internal  flow  and  the  free 
surface,  and  thus,  the  scheme  has  been  found  to  be  stable  so  long  as  the 
volume  of  water  entering  a  finite  difference  cell  within  a  time  step  is 
less  than  the  volume  of  the  cell. 

202.  Edinger  and  Buchak  utilize  the  layer-averaged  concept  of 
Leendertse  and  Simons.  The  governing  equations  that  are  solved  are  thus 
laterally  and  layer-averaged  2-D  equations  with  layer-averaged  variables 
as  the  dependent  variables.  The  equations  are  written  in  the  conserva¬ 
tive  form  with  the  Boussinesq  and  hydrostatic  approximations.  In  addi¬ 
tion,  eddy  coefficients  are  utilized  to  model  the  influence  of 
turbulence. 

203.  The  horizontal  coefficients  of  eddy  viscosity  and  eddy  dif- 
fusivity  are  assumed  to  be  constant;  whereas,  in  a  recent  development, 
the  vertical  eddy  diffusivity  and  eddy  viscosity — related  to  the  in¬ 
ternal  friction  coefficient  that  results  from  the  layer  averaging  and 
replaces  vertical  viscous  terms  as  related  to  second  derivatives — are 
allowed  to  be  dependent  upon  the  Richardson  number.  The  form  of  this 
functional  dependence  is 


Av  »  A  (1  + 
o 

e  =  e  (1  +  10R,  )-1/^2 
v  v  i 

o 

Unstable  stratification  is  modeled  by  allowing  e  to  increase  to  the 

?  v 

diffusive  stability  limit  of  Az  /2At  when  R^  <  0  . 

20l*.  As  in  other  hydrostatic  models,  the  vertical  component  of 
the  velocity  is  obtained  from  the  incompressibility  equation,  with  the 
solution  beginning  at  the  bottom  and  progressing  up  the  column  of 
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layers.  An  equation  for  the  water  surface  elevation  is  then  obtained 
by  s umning  the  layer-averaged  incompressibility  equation  over  the  water 
column.  This  equation  takes  the  form 


3(5B)  _  Y  3 (uEH ) 

3t  layers  3x 

where  t,  is  the  deviation  from  the  top  of  the  top  layer  of  fluid , 
positive  downward.  Edinger  and  Buchak  then  replace  the  time  derivative 
by  a  backward  difference  to  yield  an  implicit  solution  for  t,  .  How¬ 
ever,  the  velocities  are  unknown  at  the  n  +  1  level.  This  problem  is 
overcome  in  the  following  manner.  The  horizontal  momentum  equation 
takes  the  form 


in  which  a  forward  time  differencing  is  used  in  relation  to  all  the 
terms  comprising  F  ,  while  the  3c/3x  term  is  taken  implicitly,  i.e., 
at  the  n  +  1  time  step.  The  expression  for  (uBH)n+1  from  the  momen¬ 
tum  equation  is  then  substituted  into  the  finite  difference  form  of  the 
free  surface  equation.  The  resulting  difference  equation  then  contains 
the  unknowns  a  uridiagonal  system  results, 

which  can  be  efficiently  solved  by  che  Thomas  algorithm. 

205.  With  such  a  coupling  of  the  internal  flow  and  the  free  sur¬ 
face  computations,  the  Courant  stability  criterion  is  removed.  The  time 
step  is  now  limited  by  the  internal  flow  speed,  plus  perhaps  diffusive 
criteria,  rather  than  the  speed  of  the  surface  gravity  wave.  It  is  the 
removal  of  the  Courant  condition  that  makes  the  Edinger  and  Buchak 
(1979)  model  so  attractive  with  regard  to  long-term  simulations  of 
stratified  reservoirs. 

206.  With  the  free  surface  elevations  computed  at  the  n  +  1 
time  step,  the  horizontal  velocity  component  is  then  computed  explic¬ 
itly,  followed  by  an  explicit  computation  of  the  vertical  component. 

The  temperature  is  then  computed  from  its  transport  equation,  using 
the  new  velocities.  This,  however,  now  requires  an  implicit  solution 
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for  the  temperature  since  if  the  velocity  in  the  advective  term  is  taken 
at  the  (n+l)  level,  the  temperature  must  be  taken  at  that  level  also; 
i.e.,  terms  such  as  3(uBHT)/3x  are  taken  completely  at  the  new  time 
level.  To  avoid  having  a  2-D  implicit  computation,  which  would  require 
either  an  iterative  solution  or  perhaps  the  use  of  an  ADI  scheme,  the 
horizontal  diffusive  term  is  taken  at  the  n+l  time  level,  but  the 
vertical  diffusive  term  is  taken  at  the  old  or  n  time  level.  The  re¬ 
sulting  difference  equation  takes  a  tridiagonal  form  also  and  thus  is 
solved  in  the  same  manner  as  is  the  free  surface  equation. 

207.  Spatial  derivatives  are  replaced  by  centered  differences  in 
all  terms,  except  the  advective  terms  in  the  temperature  equation  where 
Roache's  (1972)  first  windward  differencing  is  used.  In  addition,  in 
computer  experimentation  with  the  model,  it  was  concluded  that  windward 
differencing  is  also  required  in  the  momentum  advective  terms  in  cells 
adjacent  to  forced  outlets.  This  will  be  discussed  later  in  connection 
with  application  of  the  model  to  the  GRH  flume.  Since  the  windward  dif¬ 
ferencing  is  Roache's  first  kind  in  which  simple  forward  or  backward 
differencing  is  utilized,  it  appears  the  solution  scheme  is  0(At,  Ax)  . 
As  previously  noted,  such  a  scheme  preserves  the  transportive  property 
but  not  the  conservative  property  and  in  addition  is  only  of  the  first 
order.  If  Roache's  second  upwind  differencing  had  been  employed,  the 
resulting  scheme  would  be  almost  0(At,  Ax  )  ,  and  both  the  transportive 
and  conservative  properties  would  be  preserved.  It  might  be  noted  that 
by  solving  the  temperature  equation  implicitly,  the  time  step  limit  of 
W/Az  that  can  be  more  severe  than  u/Ax  has  been  removed. 

208.  Boundary  stresses  at  the  surface  and  the  bottom  are  incor¬ 
porated  directly  into  the  layer-averaged  equations  through  the  following 
expressions : 

TtIT,lr.  =  —  p  cos  <t> 

WIND  p  Ka  a  T 


and 


t 
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where 


—  O 

C#  =  resistance  coefficient  (2.6  x  10  ) 

p  =  air  density  (1.2  kg/m  ) 

a  -1 

W  =  wind  speed  at  10-m  height  (msec  ) 

8. 

<p  =  angle  between  wind  and  reservoir  axis 

1/2 

C  =  Chezy  coefficient  (m  /sec) 

209.  With  the  use  of  the  layer-averaged  approach,  the  boundary 
stresses  are  incorporated  directly  as  terms  in  the  equations,  and  bound¬ 
ary  conditions  on  the  tangential  velocity  at  the  bottom  cannot  be  pre¬ 
scribed.  In  addition,  with  the  MAC-type  grid  employed  and  with  the 
vertical  velocity  determined  from  the  incompressibility  condition,  no 
specification  of  the  tangential  velocity  at  a  vertical  wall  is  allowed. 
Of  course,  at  all  solid  boundaries,  the  normal  component  of  velocity  is 
set  to  zero.  In  addition,  all  eddy  coefficients  are  set  to  zero  at 
solid  boundaries  to  prevent  heat  transfer  at  such  boundaries. 

210.  The  net  rate  of  surface  heat  exchange  is  expressed  by: 

h  =  -CSHE  (T  -  ET) 
n  s 

where  CSHE  and  ET  are  dependent  upon  shortwave  solar  radiation,  air 
temperature,  dew  point  temperature  and  wind  speed. 


PART  VI:  APPLICATION  OF  SELECTED  MODELS  TO  THE  GRH  FLUME 


211.  Two  of  the  3-D  models  and  three  of  the  2-D  models  have  been 
applied  to  a  bottom  density  current  problem  in  the  Generalized  Reservoir 
Hydrodynamics  (GRH)  flume  at  WES.  In  addition,  before  this  report  was 
published  Arsev  Eraslan  provided  results  from  application  of  the  auto¬ 
mat  :.c  2-D  version  of  his  general  3-D  model.*  The  3-D  models  were  that 
of  Spraggs  and  Street  (1975)  and  the  Waldrop-Tatom  (1976)  model;  while 
the  2-D  models  were  LARM,  the  TVA  model,  and  the  RMA-7  finite  element 
model.  The  two  attempts  at  a  3-D  simulation,  as  well  as  the  TVA's  2-D 
simulation,  were  made  by  the  respective  model  developers  at  the  request 
of  WES,  with  the  Waldrop-Tatom  simulation  being  made  by  Tatom  at  WES 

on  WES's  Texas  Instrument-Advanced  Scientific  Computer  (TI-ASC)  computer. 
The  simulations  with  LARM  were  conducted  by  Edinger  and  Buchak  on  the 
CDC  7600  computer  located  at  Boeing  in  Seattle,  Wash.  In  addition,  WES 
personnel  have  made  similar  computations  on  the  CYBER  176  located  at 
Kirtland  Air  Force  Base,  N.  Mex.  The  application  of  the  finite  element 
model  RMA-7  was  made  by  Bob  MacArthur  at  the  Hydrologic  Engineering 
Center  (HEC)  on  CDC  equipment  located  at  Berkeley  University  and  on  a 
Prince  550  minicomputer  located  in  Lafayette,  Calif. 

212.  The  primary  reason  for  application  of  the  models  to  the 
bottom  density  flow  problem  was  to  provide  an  assessment  of  relative 
economy  of  the  more  promising  models  and  their  ability  to  simulate  a 
real  problem  that  commonly  occurs  in  reservoirs,  whether  it  be  as  the 
result  of  a  coldwater  inflow  or  the  plunging  of  a  sediment-laden  stream. 
With  an  application  to  a  laboratory  flume,  test  conditions  could  be 
accurately  controlled  and  temperature  and  velocity  profiles  readily  ob¬ 
tained.  Although  temperature  data  are  available,  as  far  as  is  known,  a 
detailed  set  of  reservoir  field  data  including  velocities  and  results 
from  dye  tracer  tests  does  not  exist.  It  seems  reasonable  to  believe 
that  if  a  mathematical  model  can  accurately  simulate  laboratory 


*  Personal  communication,  March  1980,  Arsev  Eraslan,  Chief  Scientist, 
Hennington,  Durham,  and  Richardson,  Knoxville,  Tenn. 
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conditions,  the  expectation  oi'  reasonable  s’ i eld  abdication:;  is  justifie 
This  is  true  because  the  only  scaling  effects  in  the  mathematical  models- 
is  in  the  speci fieation  of  the  eddy  coeff icients .  Thus,  although  an 
accurate  simulation  of  a  laboratory  test  may  not  justify  a  quantitative 
confidence  in  the  ability  of  the  model  to  yield  similar  accuracy  in  the 
field,  it  does  demonstrate  qualitatively  the  model's  ability  to  simulate 
basic  flow  phenomena. 

Description  of  i'lume  and  Test  Conditions 

213.  A  photograph  of  the  JHil  flume  is  provided  in  ."inure  o.  The 
flume  is  2U.  3d  m  Ion*;  with  a  0.91-ni  *  0.9r-m  cross  section  at  the  down¬ 
stream  end.  The  cross  section  at  the  upstream  end  is  0.30  m  *  J.  30  m 


Figure  6.  The  Generalized  Reservoir  Hydrodynami  cs  (.GKii)  flu;i.e 


and  linearly  grows  in  width  over  the  first  6.10  m  to  a  cross  section 
0.30  m  deep  and  0.91  m  wide.  The  bottom  of  the  flume  is  horizontal  for 
the  first  6.10  m  and  then  drops  a  total  of  0.6l  m  linearly  over  the  fi¬ 
nal  18.29  m  of  the  flume.  Both  plan  and  side  views  are  given  in  Fig¬ 
ure  7.  The  water  in  the  flume  was  at  rest  and  homogeneous  at  the  initia¬ 
tion  of  the  test,  with  the  temperature  being  70.6°F.  Cold  water  was 
input  at  0.16  m  from  the  upstream  end  at  a  temperature  of  6 2.0°F.  A 
baffle  restricted  the  cold  water  to  enter  over  about  the  bottom  0.15  m 
of  the  cross  section.  The  inflow  rate  was  0.00063  m  /sec  with  the  out¬ 
flow  rate  at  the  downstream  end  being  the  same.  The  outflow  was  removed 
from  a  port  with  a  2.5l-cm  diameter  located  0.15  m  above  the  bottom  of 
the  flume  and  0.16  m  from  either  side.  Thus,  as  previously  discussed, 
the  2-D  laterally  averaged  models  will  not  accurately  model  the  momentum 
flux  from  the  system.  In  fact,  neither  will  a  3-D  model  unless  the  lat¬ 
eral  and  vertical  dimensions  of  a  cell  are  of  the  same  size  as  the  port. 

t 

0.91m 


a.  PLAN  VIEW 


Observed  Flow  Phenomena 


0.91m 


± 


211.  The  coldwater  input  was  dyed  for  easy  visual  observation. 
The  basic  flow  phenomena  that  developed  was  the  classical  density 
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and  turbulent  flow,  rather  than  being  completely  laminar. 

2l6.  Using  Harleman's  equation  for  the  average  velocity  of  a 
laminar  density  underflow. 


U  =  0.375  R' 


where 

h__  =  height  of  density  underflow,  ft 
Dr 

S  =  slope 

=  Reynolds  number 
Ap/p  =  0.001121 

a  value  of  U  =  0.012  m/sec  is  obtained.  If  the  equation  for  the  aver¬ 
age  velocity  of  a  turbulent  density  underflow  is  used,  i.e.. 


U  = 


Ap  hDF  5 
p  f(l  +  a) 


where  f  , 
responding 
a  =  0.1*3  ; 


the  nondimensional  friction  factor,  is  taken  as  0.003,  cor- 

1/2 

to  a  Chezy  value  of  55  m  /sec,  and  as  suggested  by  Harleman 
a  value  of  U  =  0.0U  m/sec  is  computed. 


Application  of  Three-Dimensional  Models 

217.  The  flow  in  the  flume  is  essentially  a  two-dii..~nsional  flow, 
except,  of  course,  in  the  vicinity  of  the  outlet.  However,  as  an  aid 

in  the  assessment  of  3-D  models,  an  attempt  at  applying  both  the  Spraggs 
and  Street  (1975)  and  the  nonhydrostatic  version  of  the  Waldrop- Tatom 
(1976)  models  to  the  coldwater  inflow  problem  has  been  made.  As  dis¬ 
cussed  below,  neither  of  these  attempts  was  very  successful. 

Application  of 

Spraggs  and  Street's  THERMAC 

218.  Dr.  Lynn  Spraggs  at  McGill  University  made  an  application  of 
THERMAC  with  the  computing  facilities  available  to  him  in  Montreal, 
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Canada.*  After 


uvestigation  without  success  into  the  use  of  time 


:.^s  larger  than  that  allowed  by  the  Courant  condition,  Spraggs  con¬ 
cluded  that  using  the  existing  3-D  THERMAC  model  to  simulate  the  flow 
in  the  flume  was  not  economically  feasible.  His  work  with  the  larger 
time  steps  included  different  schemes  for  temperature  acceleration, 
rigid-lid  approximations  and  differential  time-stepping.  His  estimate 
for  simulating  30  min  of  real  time  in  the  flume  is  25  to  40  hr  of  CPU 
v central  processing  unit)  time  on  a  CDC  °o0C.  Therefore ,  it  is  :.tvc  ;us 
that  numerical  schemes  that  allow  for  much  larger  time  steps  must  : 
ievisea  before  an  explicit  3-D  model  such  as  TiiERMA.i  cat.  etjts:.: 
applied  to  relatively  lone-term  reservoir  simulations. 

219.  :'■!  mss.  was  able  to  simulate  only  UjC  sec  of  the  f i  .s:.e  i  t 

lem.  The  computed  temperature  field  at  >.9  sec  is  presented  in  Pi  "u?  e 
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Figure  9.  THERMAC  Model  results  after  5**9  sec,  cooled  jet 


*  Personal  communication ,  November  1979,  Lynn  Spraggs,  Me a ill  Univer¬ 
sity,  Montreal,  Canada. 


As  can  be  seen,  the  simulation  was  conducted  with  outflow  from  the  top 
rather  than  the  bottom.  Spraggs  indicates  that  the  reason  for  the  tem¬ 
peratures  of  the  cells  at  level  11  being  greater  than  70.6°F  is  that 
the  simulation  is  unstable  at  the  surface  and  continuity  is  not  being 
conserved.  It  should  be  noted  that  in  THERMAC  the  unstable  stratifica¬ 
tion  resulting  from  the  stair-stepping  effect  at  the  bottom  is  handled 
in  a  fully  convective  manner,  since  the  complete  vertical  momentum 
equation  is  retained  and  buoyancy  effects  are  thus  convectively  modeled. 

220.  Spraggs  also  made  an  additional  simulation  with  a  heated 
bottom  inflow.  Resulting  temperatures  are  presented  in  Figure  10.  The 
simulation  shows  that  the  model  seems  to  be  performing  correctly.  How¬ 
ever,  a  much  longer  simulation  time  is  required  before  definitive  con¬ 
clusions  can  be  drawn. 


Figure  10.  THERMAC  Model  results,  heated  jet 
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Application  of 
Waldrup-Tatom  j-D  PLUME 


221.  As  previously  noted,  the  application  of  the  nonhydrustatic 
version  of  the  Waidrop-'i’atom  model  to  the  density  underflow  problem  ! u 
the  GRH  flume  was  made  by  Tatom  and  L’mith  (1979a)  on  the  i’l  ACC  computer 
located  at  WES. 

222.  The  numerical  schematization  of  the  flume  is  illustrated  in 
Figure  11.  As  shown,  the  varying  width  at  the  upstream  end  is  modeled 
with  three  regions,  each  with  a  constant  width.  The  axis  of  the  flume 
is  considered  to  be  a  plane  of  symmetry  so  that  only  half  of  the  flume 
in  the  lateral  direction  is  modeled.  As  can  be  seen  from  Figure  11,  the 
bottom  never  falls  on  a  grid  point,  and  solid  walls  are  assumed  to  lie 
halfway  between  the  last  two  rows  of  points.  Variable  grid  spacing  in 
all  three  dimensions  is  allowed  in  the  model  for  extra  flexibility. 

Very  little  documentation  of  the  code  has  been  published. 

223.  Initially,  it  was  realized  that  excessive  computing  time 
would  be  required  if  the  time  step  was  restricted  by  the  Courant  condi¬ 
tion.  With  a  lateral  spatial  dimension  of  Y.62  cm  and  a  maximum  depth 
of  0.91  m,  the  Courant  criterion  restricts  the  time  step  to  be  less  than 
approximately  0.025  sec.  Therefore,  the  initial  decision  was  made  to 
model  the  problem  using  a  rigid-lid  assumption  to  allow  for  larger  time 
steps.  Tatom  incorporated  this  by  forcing  the  water  surface  to  remain 
at  its  initial  level  and  specifying  a  derivative  boundary  condition  at 
the  surface  on  the  dynamic  pressure.  The  results  did  not  resemble  the 
density  underflow  observed  in  the  flume.  Basically,  the  ooldwater  in¬ 
flow  tended  to  spread  over  the  complete  dep>th  of  the  flume  and  no  flow 
reversal  was  computed. 

22b.  It  was  then  decided  that  the  rigid-lid  assumption  was  not 
appropriate,  and  the  derivative  pressure  boundary  condition  was  replaced 
with  a  pressure  boundary  condition  that  corresponded  to  a  free  surface. 
Actually  the  surface  was  allowed  to  be  free  only  in  the  longitudinal 
direction;  i.e.,  no  transverse  variations  were  allowed.  Applying  the 
Courant  condition  only  to  the  longitudinal  direction  gives  a  stability 
restriction  such  that  the  time  step  must  be  less  than  about  0.50  sec.  A 
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Figure  11.  Three-dimensional  grid  system 


time  step  one-tenth  of  this,  i.e.,  0.05  see,  was  then  employed.  With 
the  problem  set  up  in  this  fashion,  there  still  was  no  real  improvement 
in  the  computed  flow  field.  In  addition,  the  computing  time  was  exces¬ 
sive.  Approximately  12  to  15  hr  of  CPU  time  on  the  TI  ASC  computer  would 
have  been  required  to  simulate  30  min  of  real  time  ir.  the  flume. 

225.  Various  portions  of  the  code  were  investigated  in  an  attempt 
to  resolve  the  inability  of  the  model  to  properly  simulate  the  density 
underflow;  e.g.,  molecular  values  of  the  eddy  coefficients  corresponding 
to  laminar  flow  were  used  instead  of  the  turbulent  open  channel  coeffi¬ 
cient  model,  differencing  of  the  convective  terms  near  the  bottom  was 
changed,  and  the  pressure  boundary  condition  at  the  surface  was  modified. 
The  first  two  changes  above  made  little  or  no  difference.  When  the  pres¬ 
sure  boundary  condition  was  changed  such  that  the  dynamic  pressure  at 
the  first  row  of  grid  points  inside  the  fluid  was  set  to  zero,  some 
improvement  was  noted.  A  slight  flow  reversal  was  computed  above  the 
density  underflow.  However,  the  temperature  of  the  water  near  the  bot¬ 
tom  was  too  high  and  the  underflow  moved  much  too  slowly. 

226.  At  this  point,  Tatom  decided  again  to  invoke  the  complete 
rigid-lid  assumption  to  allow  for  a  much  larger  time  step  but  to  retain 
the  zero  dynamic  pressure  condition  at  the  surface.  Results  from  this 
run  and  a  list  of  input  parameters  are  presented  in  Appendix  A.  The 
general  conclusion  is  that  the  density  underflow  is  still  not  properly 
simulated.  As  can  be  seen  from  the  computed  results,  very  little  flow 
reversal  is  computed,  and  the  computed  flow  moves  much  more  slowly  than 
observed  in  the  flume.  Only  about  18-19  min  is  required  for  the  density 
underflow  to  traverse  the  complete  length  of  the  flume,  i.e.,  23.93  m, 
but  the  model  indicates  a  travel  distance  of  only  approximately  9.1U  m 
in  33  min.  Since  funding  provided  for  this  application  was  limited,  the 
reason  for  the  inability  of  the  model  to  properly  compute  the  density 
underflow  has  not  been  determined.  It  should  be  noted  that  in  a  mathe¬ 
matical  sense  a  Dirichlet-type  boundary  condition,  i.e.,  setting  the 
dynamic  pressure  equal  to  zero  at  the  surface,  is  not  allowed  when  impos¬ 
ing  the  rigid-lid  approximation. 

227.  As  described  in  Tatom  and  Smith  (1979b)  in  an  attempt  to 
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reduce  the  CPU  time  required  by  3-D  PLUME,  Tatom  and  Smith  recoded  por¬ 
tions  of  the  model  to  better  utilize  the  vector  features  of  the  ASC 
computer.  From  the  results  presented  in  Appendix  B,  it  can  be  seen  that 
this  effort  resulted  in  a  56-percent  reduction  in  CPU  time  over  the 
original  version  of  the  model  that  utilized  the  automated  vector  fea¬ 
tures  of  the  machine. 


Application  of  Tvo-Dimensional  Models 


Application  of 

Edinger  and  Buchak's  LARM 

228.  Because  the  Corps  funded  the  initial  development  of  LARM,  an 
early  version  of  the  basic  computer  code  was  available  for  computer  ex¬ 
perimentation  by  WES  personnel.  During  this  experimentation,  several 
general  changes  wer*  made  to  LARM.  These  centered  around  making  the 
model  more  general  in  the  specification  of  inflows  and  outflows. 

229.  During  the  application  of  LARM  to  the  GRH  flume,  as  well  as 
in  the  computer  experimentation,  it  was  observed  that  a  common  occur¬ 
rence  at  the  downstream  boundary  in  front  of  an  outlet  was  that  of  a 
flow  reversal.  Various  steps  were  taken  to  try  to  alleviate  this  prob¬ 
lem,  including  an  attempt  to  incorporate  a  momentum  correction  factor 
and  a  momentum  sink  term  under  the  assumption  that  perhaps  the  improper 
modeling  of  the  momentum  flux  through  an  outlet  was  causing  the  problem. 
In  addition,  in  an  effort  to  create  a  larger  pressure  gradient  near  the 

outlet  to  force  the  flow  in  the  proper  direction,  the  hydrostatic  pres- 

2 

sure  was  decreased  by  1/2  p  u  ,  i.e.,  the  dynamic  pressure.  None  of 
these  attempts  proved  successful. 

230.  Finally  it  was  discovered  that  the  use  of  centered  differ¬ 
ences  in  the  convective  terms  of  the  x-momentum  equation  is  unacceptable 
near  a  forced  outlet.  This  is  related  to  the  fact  that  centered  differ¬ 
ences  do  not  possess  the  transportive  property.  This  is  illustrated  by 
the  problem  below  in  which  the  initial  flow  field  is  stationary  and  a 
forced  outflow  with  velocity  Uq  is  prescribed. 
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Using  centered  spatial  derivatives  and  a  forward  time  derivative,  the 
x-velocity  at  time  t  =  At  at  (i-l/2,j)  is 


t=At  =  t=0  At 

i-l/2,J  i-l/2,j  "  AX 


r  2  2  1 

L  i.j  i-ij 


With  the  flow  field  at  t  =  0  stationary, 


U 


t=At 
i-1/2,  J 


or,  the  initial  computation  for  the  velocity  in  front  of  the  outlet 
yields  a  flow  reversal,  i.e.. 


t=At 
i-1/2, j 


As  noted  by  Spraggs  and  Street  (1975),  the  use  of  windward  differencing 

near  an  outlet  corrects  the  problem.  Therefore,  the  original  centered 

2 

difference  representation  of  the  horizontal  advective  term  3(u  BH)/3x 
has  been  replaced  with  a  one-sided  difference  near  an  outlet. 

231.  In  the  initial  application  of  LARM  to  the  GRH  flume  by  WES 
personnel,  it  was  observed  that  the  coldwater  inflow  in  essence  moved 
to  the  dam  in  the  horizontal  plane  in  which  it  entered.  The  reason  for 
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this  was  that  in  the  original  version  of  LARM,  the  eddy  coefficients 
were  not  treated  as  functions  of  the  Richardson  number  and  thus  vertical 
variations  in  density  were  not  considered.  Edinger  and  Buchak  have 
since  modified  LARM  to  allow  for  the  Richardson  number  dependence  pre¬ 
viously  presented.  Thus,  when  an  unstable  stratification  arises,  i.e., 

<  0  ,  the  vertical  eddy  viscosity  and  diffusivity  are  increased  to 
their  maximum  values  based  upon  the  diffusive  stability  criterion.  This 
procedure  forces  either  a  maximum  diffusion  upwards  or  downwards  depend¬ 
ing  upon  whether  the  density  of  the  cell  is  less  than  or  greater  than 
the  surrounding  density.  The  results  provided  by  Edinger  and  Buchak* 
(and  presented  in  Figures  16-30)  were  obtained  from  simulations  in  a 
22.87-m  flume  rather  than  the  actual  length  of  the  GRH  flume  traversed 
by  the  underflow,  i.e.,  23.93  m.  Values  of  the  various  coefficients 
and  other  input  parameters  are  presented  in  Table  2.  A  longitudinal 

Table  2 

LARM  Input  for  GRH  Flume  Application 


Parameter _  Value 


Spatial  step 

Ax  =  1.52k  m 

Layer  thickness 

H  =  0.0762  m 

Time  step 

At  =  5-0  sec 

Horizontal  viscosity 

ZT  p 

1.5  x  10  m  /sec 

Horizontal  diffusivity 

l.k  x  10  5  rn^/sec 

Vertical  viscosity  at 

1.5  x  10  m  /sec 

neutral  stability 

Vertical  diffusivity  at 

l.k  x  10  ^  m^/sec 

neutral  stability 

7n  1/2, 

70  m  /sec 

Chezy  coefficient 

Inflow 

O.OOO63  m'Vsec 

Out  flow 

O.OOO63  m  /sec 

Inflow  temperature 

62. 0°F 

*  Personal  communication,  November  1979>  J.  E.  Edinger  and  E.  M.  Buchak, 
J.  E.  Edinger  Associates,  Inc.,  Wayne,  Penn. 
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DISTANCE  OF  UNDERFLOW  FRONT  FROM  UPSTREAM.  M 


outflow  temperature. 


and  those  recorded,  tends  to  substantiate  the  conclusion  that  the  density 
underflow  is  not  being  forced  to  "hug"  the  bottom  enough.  The  stair- 
stepped  bottom  appears  to  result  in  too  much  mixing  of  the  coldwater 
underflow  with  the  warmer  water  lying  below  the  next  stair-step,  which 
results  in  higher  computed  outflow  temperatures  than  those  recorded.  The 
computed  2-D  mass  flux  field  for  60  min  after  initiation  of  the  inflow 
at  U-min  increments  is  presented  in  Figures  16-30. 

233.  It  should  be  realized  that  the  above  problem  in  the  modeling 
of  the  density  underflow  is  not  unique  to  LARM.  Any  model  that  represents 
the  bottom  boundary  in  such  a  stair-stepping  fashion  will  encounter  the 
same  problem  of  too  much  mixing  and  a  resulting  slower,  thicker,  and 
warmer  density  underflow. 

23^.  As  a  final  note,  the  results  presented  here  were  computed 
with  windward  differencing  of  the  convective  terms  throughout  the  flow 
field.  Much  smoother  computations  were  realized  than  when  centered 
differences  were  used  everywhere  except  near  the  outlet.  A  comparison  of 
the  relative  magnitude  -f  various  terms  in  the  horizontal  momentum  equa¬ 
tion  revealed  that  the  convective  terms  are  approximately  the  same  magni¬ 
tude  as  the  density  gradient  terms.  In  real  reservoirs,  convective  terms 
usually  dominate  only  in  the  backwater.  In  addition  to  the  windward 
differencing  being  employed,  the  upstream  boundary  condition  was  modified 
to  force  the  temperature  in  the  most  upstream  column  to  remain  at  the 
upstream  temperature  of  62°F,  which  resulted  in  a  slightly  faster  under¬ 
flow  current. 

Application  of  Waldrop’s  TVA  Model 

235.  Dr.  Bill  Waldrop  and  Walter  Harper  at  the  Tennessee  Valley 
Authority  in  Norris,  Tenn. ,  have  made  an  application  of  the  2-D  reser¬ 
voir  model  to  the  density  underflow  problem  in  the  GRH  flume  on  the 
computing  facilities  available  to  TVA  and  have  provided  results  to  WES.* 
Two  different  runs  were  made.  The  first  allowed  the  heavier  inflow  to 
seek  its  own  level  at  the  upstream  boundary,  but  in  the  second  run,  the 


*  Personal  communication,  November  1979,  B.  Waldrop  and  W.  Harper,  Ten¬ 
nessee  Valley  Authority,  Norris,  Tenn. 
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Computed  mass  flux  profile 


mass  flux  profile 


Figure  29.  Computed  mass  flux  profile  from  FARM  at 


l23 


coldwater  inflow  was  forced  to  enter  the  bottom  layer,  as  was  the  case 
in  the  Edinger  and  Buchak  application.  Figure  31  demonstrates  that  the 
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Figure  31.  Comparison  of  TV A  computed  and 
recorded  underflow  spec.* 

computed  underflow  moves  too  slowly.  Results  from  both  runs  Indicate  a 
computed  travel  time  in  excess  of  2h  min,  although  it  siiouLd  be  n-.de  i 
that  the  travel  time  for  the  underflow  to  traverse  the  horizontal  port  i. 
of  the  flume  agrees  quite  well  with  recorded  results.  dice  again  it 
would  seem  these  results  tend  to  substantiate  tiie  previous,  statements 
made  concerning  the  stair-stepping  effect  of  the  bottom.  This  effect  is 
further  indicated  from  the  plot  of  computed  versus  recorded  outflow  tem¬ 
peratures  presented  in  Figure  32.  Computed  velocity  fields  and  isothern 


Figure  32.  Comparison  of  TVA  computed  and 
recorded  outflow  temperature 

from  the  first  application  are  presented  in  Figures  33-37  at  times  of 
6,  12,  18,  2h,  and  30  min  after  initiation  of  the  inflow.  Similar  plots 
from  the  second  application  are  presented  at  10,  20,  and  30  min  in  Fig¬ 
ures  38-1*0.  The  only  results  provided  for  a  direct  comparison  of  computed 
and  recorded  velocities  at  a  particular  location  are  presented  in  Fig¬ 
ure  1*1.  There  a  comparison  of  the  computed  velocities  at  10. 67  m  from 
the  upstream  end  at  10  min  after  initiation  is  made  with  recorded  veloc¬ 
ities  at  11.1*3  m  from  the  upstream  end  at  11  min  after  initiation. 

236.  As  previously  noted,  the  Waldrop  model  is  an  explicit  FDM 
and  thus  the  time  step  at  which  computations  are  male  is  restricted  by 
the  speed  of  the  surface  gravity  wave.  For  thio  pplication,  the  maximum 
allowable  time  step  is  computed  to  be  about  0.30  sec.  Waldrop  indicates 
that  a  time  step  of  0.30  sec  was  actually  used,  which  resulted  in  1*6  sec 
of  CPU  time  on  a  CDC  7600  computer  for  6000  time  steps. 

237.  As  a  final  note,  Waldrop  has  indicated  that  he  also  has  en¬ 
countered  flow  reversals  in  front  of  forced  outlets.  However,  rather 


125 


20.0  « 


r 


Figure  hi.  Comparison  of  TVA  computed  and 
recorded  horizontal  velocities 


than  using  windward  differencing  near  the  outlet,  he  reduces  the  pres¬ 
sure  by  subtracting  the  dynamic  pressure  l/2plf  where  U_  is  the 

G  G 

outlet  velocity,  to  force  the  flow  in  the  proper  direction.  This  is 
easy  to  implement,  since  his  basic  computational  ceil,  as  is  illustrated 
below,  has  the  velocity  components  defined  at  the  center  of  the  cell 
with  pressures  defined  on  the  vertical  faces. 


Application  of 

Norton-King-Orlob  FEM — RMA-7 

238.  The  application  of  the  Norton-King-Orlob  FEM  to  the  density 
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underflow  in  the  GRH  flume  was  made  by  Bob  MaeArthur  of  the  Hydrologic 
Engineering  Center  (HEC)  in  Davis,  Calif.*  As  in  the  previous  model 
applications,  the  initial  conditions  consisted  of  zero  flows  (zero  veloc¬ 
ity  throughout)  and  isothermal  water  temperatures  throughout  at  T0.6°F. 
After  time  zero,  a  constant  coldwater  inflow  of  O.OOO63  m  /sec  at 
62. 0°F  was  imposed  entering  near  the  bottom  of  the  flume  as  an  upstream 
boundary  condition.  MaeArthur  indicated  that  a  zero  pressure,  free  dis¬ 
charge  boundary  condition  was  prescribed  at  the  outlet  so  the  inflow 
rate  would  equal  the  outflow  and  the  free  water  surface  would  remain 
horizontal.  Values  of  the  eddy  coefficients  used  are  presented  in 
Table  3. 


Table  3 

Values  of  the  Turbulent  Exchange  Coefficients 


Used  for  the  GRH  Flume  Applications 


Turbulent  Exchange  Coefficients 


(Eddy  Viscosity 


Value 


£ 

XX 

e 

xy 

e 

yx 


10  lb-sec /ft2  (0.48  m2/sec) 

0.004  lb-sec/ft2  (1.9  *  10  ^  m2/sec) 
10  lb-sec/ft2  (0.48  m2/sec) 

0.01  lb-sec/ft2  (4.8  x  10  ^  m2/sec) 


Turbulent  Diffusion  Coefficients 


(Eddy  Diffusivit^ 


D 

x 

D 

y 


_ Value _ 

2.0  ft2/ sec  (0.19  m2/ sec ) 

0.01  ft2/sec  (9-3  x  10-1"  m2/sec) 


239.  Using  a  time  step  of  3  min,  a  total  time  of  18  min  was 
simulated.  Results  furnished  by  MaeArthur  for  the  first  three  time  steps 
are  presented  in  Figures  42,  43,  and  44.  The  finite  element  network 
used  to  simulate  these  results  is  shown  in  Figure  45.  Although  compara¬ 
tive  plots  are  not  presented,  MaeArthur  stated  that  travel  times  for  the 


*  Personal  communication,  November  1979,  R.  C.  MaeArthur,  Hydrologic 
Engineering  Center,  Davis,  Calif. 
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Figure  1+2.  Velocities  computed  by  RMA-T  at  T  - 
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Figure  1+3.  Velocities  computed  by  RMA-T  at  T 
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Figure  1*1+.  Velocities  computed  by  RMA-7  at  T  =  9  min 


RMA7  APPLICATION  TO  THE  WES  GRH  FLUME  DATA  JULY  1979  BY  RC  MACARTHUR 


Figure  U5.  Finite  element  network  for  GRH  flume  application 
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coldwater  plume  to  reach  the  outlet  were  usually  about  18  to  20  min  from 
the  time  of  initial  inflow.  However,  outlet  temperatures  were  several 
degrees  (°F)  greater  than  the  temperatures  measured  in  the  physical  model 
tests.  In  addition,  the  thickness  of  the  plume  and  temperatures  within 
the  bottom  flow  plume  were  greater  than  observed.  This  indicates  more 
vertical  mixing  is  occurring  in  the  numerical  simulation  than  was  occur¬ 
ring  in  the  GRH  flume,  which  MacArthur  attributes  to  the  choice  of  the 
vertical  mixing  coefficients. 

2h0.  An  inspection  o^  the  two-dimensional  velocity  field  pre¬ 
sented  after  times  of  3,  6,  and  9  min  in  Figures  42,  1* 3 ,  and  44,  respec¬ 
tively,  reveals  that  the  classical  density  underflow  phenomenon  has  not 
fully  developed  in  the  computed  results.  At  the  ponnt  where  the  flume 
bottom  begins  sloping,  the  flow  seems  to  he  projected  across  in  a  hori¬ 
zontal  plane.  One  reason  for  this  may  be  the  large  values  assumed  for 
some  of  the  eddy  coefficients,  e.g.,  the  diagonal  component  in  the  x- 

direction  e  =  0.48  m  /sec  .  MacArthur  has  indicated  that  the  stability 
XA 

of  the  model  is  extremely  sensitive  to  the  values  used  for  these  coeffi¬ 
cients.  Therefore,  such  large  values  were  required  to  obtain  stable 
solutions . 

241.  The  computer  time  on  a  CDC  76OO  to  simulate  the  results  pre¬ 
sented  here  required  h2  sec  of  execution  time  to  compute  18  min  of  flow 
time.  This  compares  with  the  approximately  5  sec  required  by  LARM  to 
compute  30  min  of  flow  time  and  46  sec  required  by  the  Waldrop  explicit 
model  to  also  simulate  30  min  of  flow  time. 

242.  After  these  initial  runs,  MacArthur  made  some  comparative 
runs,  using  the  GRH  flume  geometry,  for  homogeneous  flow  conditions  and 
additional  thermally  stratified  flow  conditions.  In  each  case,  flows  of 
O.OOO63  m  /sec  were  introduced  with  a  linear  velocity  distribution  in 
the  bottom  element  at  the  upstream  end  of  the  flume.  Figure  45  presents 
the  finite  element  network  used  for  the  simulations.  The  homogeneous 
case  was  run  isothermally  at  a  temperature  of  50.5°F,  while  the  nonhomo- 
geneous  case  was  started  with  an  initial  temperature  of  50.5°F  throughout 
and  an  inflowing  water  temperature  of  4l.0°F. 

243.  Velocity  distributions  produced  by  these  comparative  runs 
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are  presented  in  Figure  1+6  for  three  different  time  steps  (after  3,  9, 
and  18  min)  at  the  two  sampling  stations  indicated  on  Figure  1+5-  The 
effects  of  flow  stratification  are  quite  evident  and  appear  to  be  quali¬ 
tatively  reasonable. 

Application  of 

Eraslan's  Discrete  Element  Model 

21+1+ .  After  the  initial  writing  of  this  report,  Eraslan  provided 
results  from  applications  of  his  model  called  FLOWER.*  FLOWER  is  a 
computer  code  for  simulating  fast-transient  three-dimensional  coupled 
hydrodynamic,  thermal  and  salinity  conditions  in  the  intake  and  dis¬ 
charge  zones  of  power  plants  operating  on  rivers,  lakes,  estuaries  and 
coastal  regions.  The  general  3-D  model  contains  an  automatic  2-D  later¬ 
ally  averaged  version,  which  was  used  in  the  GRH  flume  simulations. 

2l+5-  Eraslan  indicates  that  the  turbulent  transport  model  of 
FLOWER  is  completely  closed;  i.e.,  it  utilizes  the  same  turbulent  (and 
laminar)  transport  model  for  all  time  and  spatial  scales  in  applications 
to  vastly  different  problems,  including  the  scales  of  physical  models 
as  well  as  the  scales  of  prototype  conditions.  Therefore,  the  user 
never  specifies  any  friction  or  turbulent  diffusion  coefficients . 

246.  Two  separate  simulations  were  made  with  the  flume  discretized 
as  shown  in  Figure  47.  One  was  a  coldwater  inflow,  while  the  other  was 
a  hotwater  input.  The  coldwater  inflow  simulation  was  the  same  as  pre¬ 
viously  discussed,  with  the  exception  that  the  outflow  was  0.00109 
m  /sec  and  the  inflow  temperature  was  54.l4°F.  Therefore,  the  water 
surface  dropped  slightly  during  the  simulation.  Figures  48-55  present 
"snap  shots"  at  200-sec  intervals  of  the  velocity  field  with  no  exagger¬ 
ation  of  the  vertical  component  for  1600  sec  after  initiation  of  the 
inflow.  From  an  inspection  of  Figure  52,  it  can  be  seen  that  the  com¬ 
puted  travel  time  required  to  traverse  the  flume  is  16-17  min,  which  com¬ 
pares  reasonably  well  with  a  recorded  time  of  approximately  15  min  for 
these  input  conditions.  Eraslan  indicates  that  if  the  inflow  had  been 


* 


Personal  communication,  March  1980,  Arsev  Eraslan,  Chief  Scientist, 
Hennington,  Durham,  and  Richardson,  Knoxville,  Tenn. 
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simulated  velocity  distributions 


METRES 


Figure  U7.  Schematization  of  GRH  flume  in 
Eraslan's  application 

1.01  M/SEC  TIME  *  200.0  SEC 
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Figure  48.  Velocities  computed  by  Eraslan's  Model  at 
T  =  200  sec,  coldwater  inflow 
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Figure  1+9.  Velocities  computed  by  Eraslan's  Model  at 
T  =  1+00  sec ,  coldwater  inflow 
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Figure  c.0.  Velocities  computed  by  Eraslan's  Model  at 
T  =  600  sec,  coldwater  inflow 
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Figure  51.  Velocities  computed  "by  Eraslan's  Model  at 
T  =  800  sec ,  coldwater  inflow 
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Figure  52.  Velocities  computed  by  Eraslan's  Model  at 
T  =  1000  sec,  coldwater  inflow 
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Figure  55<  Velocities  computed  by  Eraslan's  Model  at 
T  =  1600  sec,  coldwater  inflow 

specified  over  only  the  bottom  half  of  the  cross  section  rather  than 
the  complete  section,  a  faster  and  less  thick  underflow  current  would 
have  resulted. 

2i*7.  As  noted  above,  an  additional  simulation  was  made  in  which 
warm  water  at  70.6°F  was  input  uniformly  over  the  upstream  end  with  the 
water  in  the  flume  initially  being  stationary  and  homogeneous  at  a 
temperature  of  5^.1^°F.  Figures  56-65  present  "snap  shots"  at  200-sec 
intervals  of  the  resulting  2-D  flow  field  for  2000  sec  after  initiation 
of  the  inflow. 

2U8.  Since  FLOWER  is  an  explicit  model,  the  time  step  is  re¬ 
stricted  by  the  gravity  wave  stability  criterion  based  upon  the  deepest 
part  of  the  flume.  The  results  presented  were  obtained  by  Eraslan  from 
running  FLOWER  on  a  PDP-10  computer. 
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Figure  56.  Velocities  computed  by  Eraslan's  Model 
at  T  =  200  sec,  warmwater  inflow 
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Figure  57.  Velocities  computed  by  Eraslan's  Model 
at  T  =  U00  sec,  varmvater  inflow 
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igure  58.  Velocities  computed  by  Eraslan's  model  at  T  =  600  sec, 

warmvater  inflow 
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Figure  59.  Velocities  computed  by  Eraslan’s  model  at  T  =  800  sec, 

warmwater  inflow 
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Figure  60.  Velocities  computed  by  Eraslan's  Model 
at  T  =  1000  sec,  warmwater  inflow 
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Figure  6l .  Velocities  computed  by  Fr.arlan '  s  Model 
at  T  =  1200  sec,  warmwater  inflow 
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Figure  62.  Velocities  computed  by  Eraslan' s  Model 
at  T  =  lUOO  sec,  warmvater  inflow 
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Figure  63.  Velocities  computed  by  Eraslan' s  Mode; 
at  T  =  l600  sec,  warmvater  inflow 
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Figure  C4.  Velocities  computed  by  Eraslan's  Med 
at  T  =  lSOG  sec,  wamwater  inflov 
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65.  Velocities  computed  by  Eraslan's 
at  T  =  2000  sec,  warmwater  inflow 
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PART  VII:  CONCLUSIONS  AND  RECOMMENDATIONS 


2l*9-  Many  different  types  of  numerical  hydrodynamic  models  exist. 
These  range  from  steady  to  unsteady  models  with  the  physical  problem  area 
represented  by  one,  two,  or  three  spatial  dimensions.  In  addition,  some 
models  consider  the  effect  of  temperature  and/or  salinity  on  the  density 
of  the  water;  whereas,  others  treat  the  water  body  as  being  homogeneous. 
Most  numerical  hydrodynamic  models  invoke  the  Boussinesq  approximation 
as  well  as  the  hydrostatic  pressure  assumption;  however,  there  are 
models  that  do  neither  and  are  thus  able  to  convectively  model  buoyancy 
effects.  Some  models  allow  for  the  movement  of  a  free  surface  and  its 
subsequent  effect  on  the  internal  flow;  whereas,  others  impose  a  mathe¬ 
matical  rigid-lid  approximation  to  enable  larger  time  steps  to  be 
employed  in  the  numerical  solution  technique.  A  vast  majority  of  the 
hydrodynamic  models  employ  the  finite  difference  method  to  develop  nu¬ 
merical  solutions,  although  there  are  existing  models  that  employ  the 
finite  element  method  for  the  spatial  integration  of  the  governing 
fluid  dynamic  equations.  The  vast  majority  of  numerical  hydrodynamic 
models  handle  the  exchange  of  energy  from  the  large-scale  circulation 
patterns  to  the  small-scale  unresolvable  eddies  through  the  use  of  eddy 
viscosity  and  diffusivity  coefficients.  However,  there  are  substantial 
differences  in  the  expressions  used  to  relate  these  eddy  coefficients 
to  properties  of  the  mean  flow  field. 

250.  One-dimensional  models  are  often  applied  to  reservoirs  where 
the  principal  variation  cf  flow  characteristics  is  in  the  vertical 
direction.  The  primary  advantage  of  such  models  is  their  ability  to 
resolve  long-term  or  seasonal  temperature  profiles  economically.  Such 
models,  however,  are  not  applicable  for  predicting  multidimensional  flow 
fields  within  stratified  reservoirs  for  quality  predictions.  Therefore, 
only  two-  and  three-dimensional  models  have  been  investigated  in  this 
study. 

Conclusions  on  Two-Dimensional  Modeling 
In  order  for  a  numerical  hydrodynamic  model  to  be  applicable 
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251. 


r 


it  must  first 


to  the  prediction  of  flov  fields  in  stratified  reservoirs, 
of  all  be  at  least  a  two-dimensional  ( verti  cal-lorigi  t udinal )  model  and 
preferably  one  that  is  laterally  averaged  to  account  for  width  chan rr-:: 
along  the  axis  of  the  reservoir  as  well  as  with  depth.  The  r.odel  rust 
be  dynamic,  i.e.,  time-dependent,  and  must  be  a  heat-conducting  r.odel 
that  can  handle  unstable  stratifications.  In  other  words,  surface  hea4 
exchange  and  a  subsequent  modeling  of  the  temperature  field  and  its 
coupling  with  the  flow  field  through  its  influence  on  the  wafer  density 
must  be  handled.  In  addition,  since  the  model  will  be  applied  ever 
natural  stratification  cycles,  during  which  significant  flooding  car. 
occur,  a  free  surface  must  be  allowed  as  cpjposed  tc  the  rigid-oid  ar;  rex- 
imation.  These  are  necessary  criteria.  Considerations  of  accuracy  and 
economy  must  naturally  be  taken  into  account  also  when  selecting  a  model 
for  widespread  use  throughout  the  Corps. 

252.  Of  the  various  two-dimensional  models  investigated,  six  f 

models  come  close  to  meeting  the  required  criteria  outlined  above.  The.ce 

are  the  models  of  Edinger  and  Buchak;  Waldrop;  Thompson;  Norton,  King, 
and.  Orlob ;  Roberts  and  Street;  and  Slotta  et  al.'s  NUMAC.  The  criteria 
satisfied  by  these  models  are  summarized  in  Table  It .  In  addition, 
although  an  in-depth  investigation  has  not  been  made  due  tc  a  lack  rf 
published  material  as  well  as  publication  deadlines  tc  be  met,  resultr 
from  Eraslan's  2-D  simulations  imply  that  it  also  meets  these  criteria. 

253.  Thompson's  laterally  averaged  model  is  being  developed 

primarily  for  near  field  selective  withdrawal  studies.  The  governing 
equations  are  solved  implicitly  using  an  iterative  technique.  Thus, 
although  the  model  will  be  a  completely  general,  fully  convective  model 
that  will  accurately  handle  general  boundaries  through  the  use  of 
bound -1  ^ed  coordinates,  the  computing  costs  for  long-term  simula¬ 

tions  will  probably  prohibit  its  use  over  natural  stratification  cycles 
in  reservoirs. 

25^4.  The  Roberts  and  Street  model  assumes  a  hydrostatic  pressure 
but  does  handle  unstable  densities  by  allowing  for  a  large  diffusion  of 
heat  within  the  unstable  water  column.  Of  the  six  2-D  models  that 
satisfy  the  necessary  criteria,  this  model  and  NUMAC  are  the  only 
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Comparisons  of  Different  Models 


rH  G 

a  so  o 

O  G  *H 
•H  Q  P 

a  p  ^ 

O  H 

P  P  3 
O  o  s 
a  -h 

w  co 


P 

0 

1 — 1 

P 

P 

P 

P 

a 

<D 

o 

<D 

cd 

o 

&> 

1 

<L> 

P< 

P 

1 — 1 

p 

P 

tiO 

P 

P 

o 

<u 

<L> 

oJ 

P 

c 

P 

cd 

G 

W 

EH 

CO 

*«H 

CO 

rd 

PQ 

0 

P 

1 

1 

c 

fcO 

G 

1 

i 

1 

S 

p 

Pv 

w 

o 

P 

i 

O 

10 

Ph 

P 

w 

<D 

O 

bO 

p 

P 

c 

to 

p 

o 

(V 

G 

TJ 

P 

S0 

rH 

P 

o 

o 

P< 

p 

p 

bD 

2 

G 

Tj 

cd 

•H 

Cd 

p 

a 

<u 

T3 

C 

a 

(U 

o 

rH 

P 

g 

p 

p 

p 

o 

p 

i — i 

•rH 

•rH 

<D 

•H 

Cti 

Ph 

i — 1 

3 

0 

rC 

o 

cd 

rJ 

CO 

►P 

►P 

CO 

PQ 

a 

Eh 

« 

W 

pure  2-D  models.  The  other  four  solve  the  laterally  averaged  equations. 
The  major  disadvantage  of  the  Roberts  and  Street  model  lies  in  its  use 
of  a  numerical  solution  technique  that  restricts  the  time  step  to  be 
smaller  than  the  time  required  for  the  surface  gravity  wave  to  traverse 
a  computational  cell.  Such  a  restriction  can  result  in  excessive  com¬ 
puter  costs  for  applications  extending  over  several  months. 

255*  The  Norton,  King,  and  Orlob  model  is  similar  to  the  Thompson 
model  in  that  the  complete  vertical  momentum  equation  is  solved. 

However,  the  Norton,  King,  and  Orlob  model  uses  the  finite  element 
method  to  perform  the  spatial  integration  of  the  governing  equations. 
With  the  use  of  the  finite  element  method,  boundary  geometry  can  be 
accurately  handled  but  computing  costs  may  become  excessive  for  long¬ 
term  simulations.  Another  disadvantage  is  that  although  a  free  surface 
is  allowed,  modifications  would  probably  be  required  to  allow  for  large 
fluctuations  that  might  occur  over  a  stratification  cycle.  Wi  h  the 
complicated  coding  of  finite  element  models,  it  appears  that  £  lificant 
modifications  can  often  become  major  tasks. 

256.  The  NUMAC  model  is  based  upon  the  MAC  work  1  '  Welch  et  al , 
and  as  such,  like  the  Thompson  and  Norton,  King,  and  Orlob  model;,  is 
a  completely  convective  model.  Once  again,  however,  the  computing  costs 
for  long-term  simulations  would  be  excessive.  Not  only  is  a  two- 
dimensional  Poisson  equation  solved,  but  the  basic  computations  utilise 
an  explicit  solution  technique  with  the  maximum  time  step  restricted 
by  the  speed  of  the  surface  gravity  wave. 

257*  The  2-D  Waldrop  model  appears  to  be  a  well  developed 
laterally  averaged  hydrostatic  model  that  can  be  directly  applied  in 
its  present  form  to  predict  stratified  reservoir  hydrodynamics.  The 
manner  in  which  the  bottom  boundary  condition  on  the  velocity  is 
prescribed  would  seem  to  allow  for  a  more  accurate  modeling  of  flow 
near  the  bottom.  However,  the  bottom  is  still  in  essence  represented 
in  a  stair-step  fashion,  which  results  in  excessive  mixing  of  density 
underflows.  This  can  be  seen  from  the  results  of  the  application  to 
the  GRH  flume.  The  major  disadvantage  of  the  Waldrop  model  (and 
Eraslan’s  discrete  element  model)  is  the  gravity  wave  restriction  on 
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the  time  step  as  a  result  of  the  explicit  finite  difference  scheme  em¬ 
ployed.  As  Waldrop  has  noted,  if  the  boundary  conditions  are  varying 
rapidly  enough  to  require  input  at  time  intervals  on  the  order  of  the 
maximum  time  step  allowable  by  the  Courant  condition,  explicit  models 
can  often  be  shown  to  be  more  economical  than  implicit  ones  due  to  their 
less  complicated  coding.  However,  if  an  extremely  general  model  is 
desired  for  use  in  long-term  reservoir  simulations  during  which 
boundary  input  may  or  may  not  be  rapidly  varying,  it  appears  difficult 
to  justify  the  selection  of  a  model  with  the  time  step  restricted  by 
the  Courant  condition. 

253.  The  Hunger  and  Buchak  model  (LARM)  is  a  laterally  averaged, 
hydrostatic  model  that  employs  a  unique  method  for  removing  the  Courant 
condition  as  a  stability  criterion.  This  is  accomplished  through  a 
coupling  of  the  water  surface  computations  and  the  internal  flow  such 
that  the  water  surface  is,  comf  uted  implicitly,  while  the  internal  com¬ 
putations  are  5  er formed  explicitly.  Unstable  stratifications  arc  in¬ 
directly  handled  ty  forcing  the  maximum  diffusion  allowed  by  the 
stability  criterion  into  adjacent  cells.  Results  from  applications  of 
both  the  Edinger  and  Buchak  and  the  Waldrop  2-D  models  to  the  GRH  flume 
are  encouraging.  in  addition,  the  results  from  the  2-D  version  of 
Eraslan's  3-D  code  agreed  quite  well  with  the  flow  phenomena  observed 
in  the  flume  for  his  input  conditions. 

259-  There  are  several  areas  of  the  Edinger  and  Buchak  model 
that  should  be  investigated  for  possible  further  development.  With  its 
modular  programming,  significant  modification  of  the  model  should  not 
be  unduly  difficult.  These  areas  are  discussed  later.  Because  the 
Edinger  and  Buchak  model  satisfies  the  necessary  criteria — namely, 
time-dependent,  free  surface,  2-D  laterally  averaged,  variable  density 
and  heat-conducting — and  allows  for  unstable  stratification  and  the 
solution  technique  allows  for  economical  long-term  simulations,  it  is 
the  most  logical  2-D  model  to  select  for  further  development  to  provide 
the  Corps  with  an  accurate  and  economical  predictive  capability  in  the 
area  of  reservoir  hydrodynamics. 
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Conclusions  on  Three-Dimensional  Modeling 


260.  The  state  of  the  art  is  such  that  it  does  not  appear  any  of 
the  three-dimensional  models  investigated  can  be  economically  applied 
for  long-term  reservoir  simulations.  However,  since  most  reservoirs 
actually  exhibit  a  three-dimensional  nature,  undoubtedly  the  need  within 
the  Corps  for  a  three-dimensional  predictive  capability  will  increase 
over  the  next  few  years.  To  satisfy  this  need  in  a  practical  sense, 
new  solution  techniques  as  well  as  increased  computing  power  must  be 
realized.  In  addition,  one  should  consider  making  the  hydrostatic  pres¬ 
sure  assumption  to  remove  the  computing  cost  of  solving  for  a  nonhydro¬ 
static  pressure. 

261.  Neither  of  the  3-D  models  applied  to  the  GRH  flume  yielded 
very  encouraging  results.  Spraggs  was  not  able  to  simulate  more  than 
600  sec  with  THERMAC  due  to  the  extremely  long  computing  times  required. 
The  nonhydrostatic  version  of  the  Waldrop-Tatom  model  was  run  with  the 
free  surface  frozen,  which  allowed  a  large  time  step  to  be  used.  How¬ 
ever,  the  density  underflow  was  not  properly  computed.  After  33  min, 
the  model  computed  a  travel  of  only  approximately  10. 06  m;  whereas,  only 
18  min  was  required  in  actuality  for  the  underflow  to  traverse  the 
complete  length  of  the  flume  (23.93  m).  Tatom  feels  that  the  problem 

is  related  in  some  manner  to  the  dynajnic  pressure  computations.  There¬ 
fore,  if  time  and  funds  had  permitted,  it  would  have  been  interesting 
to  apply  the  hydrostatic  version  to  the  same  density  underflow  problem. 

Recommendations  for  Two-Dimensional  Modeling 

262.  As  noted  above,  it  is  believed  that  the  Edinger  and  Buchak 
2-D  laterally  averaged  model  offers  the  most  promise  in  the  area  of 
multidimensional  stratified  reservoir  hydrodynamic  modeling  in  the  near 
future.  However,  additional  developmental  work  and  modifications  are 
needed  to  make  the  model  more  flexible  and  accurate;  therefore,  it  is 
recommended  that  the  following  items  be  investigated  during  the  next  year 
for  possible  further  development  and  incorporation  into  HARM: 


a.  Thi  application  to  the  density  underflow  in  the  GRH 
flume  demonstrates  the  mixing  effect  of  a  stair- 
stepped  bottom.  It  is  believed  that  a  transformation 
of  the  vertical  coordinate,  as  is  performed  in  Lick's 
model,  offers  one  solution  to  this  problem.  Particular 
attention  should  be  directed  toward  the  neglection  of 
the  cross-derivative  terms  resulting  from  the  non- 
orthogonal  transformation. 

b_.  LARM  presently  allows  for  the  vertical  eddy  coefficients 
to  be  functions  of  the  Richardson  number;  however,  tne 
horizontal  coefficients  are  assumed  to  be  constant.  It 
is  recommended  that  an  eddy  coefficient  model  similar 
to  that  of  Spraggs  or  perhaps  the  simpler  model  employed 
by  Waldrop  and  Harper  (or  perhaps  Eraslan's  closed 
turbulence  model)  be  incorporated  into  LARM.  This 
should  be  relatively  easy  to  accomplish,  since  the  com¬ 
puter  code  was  initially  programmed  with  such  an 
addition  in  mind. 

c_.  LARM  presently  employs  either  windward  or  centered  dif¬ 
ferences  to  represent  the  advective  terms  in  both  the 
momentum  and  the  temperature  transport  equations  along 
with  a  forward  time  difference.  Such  a  first  order 
transport  scheme  is  adequate  for  continuous  distribu¬ 
tions.  However,  if  instead  of  an  essentially  continuous 
distribution,  a  slug  of  some  quality  constituent  is  to 
be  traced  through  the  reservoir,  large  errors  can  result 
from  the  use  of  such  first  order  schemes.  Therefore, 
it  is  recommended  that  higher  order  transport  schemes 
such  as  those  employed  by  DHI  or  perhaps  the  2-point 
scheme  of  Holly  and  Preissman  be  investigated  for  use 
in  the  modeling  of  quality  constituents,  rather  than 
the  scheme  LARM  presently  uses  for  temperature. 

cl.  It  would  seem  that  Waldrop's  method  of  setting  the 
bottom  boundary  condition  on  velocity  by  matching  a 
logarithmic  profile  is  quite  realistic.  It  is  recom¬ 
mended  that  the  use  of  such  a  boundary  condition  in 
LARM  be  investigated.  The  layer-averaged  approach 
taken  in  LARM  may  make  this  difficult. 

£.  The  horizontal  grid  spacing  in  LARM  is  constant.  A 

variable  grid  spacing  would  be  useful  to  provide  greater 
flexibility  in  the  resolution  of  a  quality  constituent 
in  a  particular  area.  The  difficulty  in  allowing  this 
and  the  subsequent  errors  that  might  occur  should  be 
determined.  Obviously,  the  linear  averaging  now 
employed  to  provide  values  of  variables  at  points  where 
they  are  not  defined  would  have  to  be  changed  to  reflect 
a  weighted  average.  Also,  as  discussed  by  Brown  and 
Pandolfo  (1979),  a  nonuniform  spatial  grid  can  influence 
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the  numerical  stability  of  the  transport  ••  ;-*ati- 
Recall  also  the  deterioration  of  ti.e  tr-ii.ou*  ion  error- 
associated  witli  variable*  ,-ri  it. 

LARM  models  a  single  reach,  i.e.,  t:.--  r--u  ..  .win,- 

major  axis  of  the  reservoir.  iesi ratio  fiexi i  i.  ify 
is  tiie  ability  to  simultaneously  model  :  ay 

of  tiie  reservoir  as  well  as  its.  major  arms  •  -r  tritu- 
taries,  i.e.  ,  to  allow  for  a  tree  structure  r  tranche-: 
network.  The  possibility  of  making  idVRM  a  multi,'  unction 
model  should  be  considered. 

The  applicability  of  LARM  at  the  downstream  •  -nu  of  a 
reservoir  near  the  dam  is  questionable  due  to  the  hydro¬ 
static  pressure  assumption  as  well  as  tiie  Inaccurate 
modeling  of  tiie  momentum  flux  through  tin-  outlet.  The 
hydrostatic  pressure  assumption  is  j  r-:  tally  as  septal  , 
except  in  tiie  inrnted  i  ate  vicinity  of  tiie  outlet;  whereas , 
tiie  inaccurate  modeling  of  the  momentum  flux  us  a 
result  of  the  model  being  less  than  three— 1  i mens iona.. 
might  extend  for  several  - i - ui.  widths.  The  ise  >f  a 
momentum  correct  i on  fact  >r  to  bettor  model  the  momentum 
flux  near  tiie  dam  should  be  investigated.  In  addition, 
the  coupling  with  1  .Ail’d  >••!*  a  model  that  is  more  appli- 
cf.it  le  near  tin  vicinity  »f  an  outlet,  e.g.  ,  Ronari  and 
Trace  (igt'd),  to  j  rov'DR-  t.he  downs. tream  boun  iary  oon- 
dition  for  LARM  ::T,-iit  be  c. ens i  i.-r- -i . 

LARM  has  been  leveloped  s  ion  taut  tiie  user  ins  to 
"hard  wire,"  i.e.,  uys. ;  cal  ]y  ciruige  code  statements , 
tiie  model  for  each  a;  plication,  ai  thougn  'lie  vers,  i  on 
being  run  by  Wi-iT  pers  ir.no  L  has.  been  modified,  to  pro¬ 
vide  a  slightly  more  general  model .  it  is  not  believed 
that  the  current  approach  of  forcing  a  user  to  change 
code  statement,  for  each  np:  locution  is  acce; table  in 
a  model  to  be  made  available  to  all  Corps  division  arid 
District  offices.  Therefore,  '.ARM  siiuud  i  be  made 
sufficiently  general  so  that  it.  can  be  a|ilied  to  a 
wide  range  of  problems  strictly  through  a  change  of 
input  data  only. 

In  the  modeling  of  real  reservoirs,  no  iften  encounters, 
side  embayments  that  contribute  es.sent  i  al  ly  nothing  t< 
tiie  momentum  of  tiie  flow  field,  but  must  still  be 
accounted  for  as  storage  areas,  in  the  conservation  of 
mass..  More  accurate  computations,  would  be  realised  i  f 
LARM  allowed  for  tiie  spec  i  f  i  cat  i  on  of  two  wi  itiis.  no 
for  tiie  total  width  that  is.  currently  input,  which 
would  continue  to  lie  used  iri  tiie  continuity  equation. 

The  other  would  be  a  width  correspond ing  to  the  actual 
flow  area  for  use  in  tiie  momentum  equation. 


j_.  In  the  use  of  numerical  hydrodynamic  models,  the  ini¬ 
tial  state  of  the  system  as  well  as  time-varying 
boundary  conditions  must  be  prescribed  as  input.  The 
nature  of  the  hyperbolic  equations  being  solved  is  such 
that  after  a  sufficient  length  of  time,  the  effect  of 
the  initial  conditions  becomes  negligible.  The  time 
simulated  to  remove  initial  condition  effects  is  com¬ 
monly  referred  to  as  the  "start-up  time."  One  way  to 
handle  this  problem  is  to  compute  a  steady  state  before 
imposing  the  time  variation  of  the  boundary  conditions. 
LARM  presently  can  only  compute  a  steady  state  as  the 
asymptotic  convergence  of  a  time-varying  solution  com¬ 
puted  by  holding  the  boundary  conditions  constant.  The 
possibility  of  incorporating  into  LARM  the  capability  of 
solving  the  steady-state  equations,  as  allowed  by  the 
Norton,  King,  and  Orlob  FEM,  should  be  investigated. 

After  the  above  modifications  are  made,  in  particular 
Lick's  transformation  to  allow  for  better  representation 
of  the  bottom,  LARM  should  again  be  applied  to  the 
density  underflow  problem  in  the  GRH  flume.  In  addition, 
hopefully  a  good  set  of  field  data  will  be  available  by 
then  through  EWQOS.  Assuming  an  eddy  coefficient  model 
similar  to  that  of  Spraggs  and  Street  (1975)  has  been 
incorporated,  these  two  applications  should  provide  in¬ 
formation  on  whether  a  single  scaling  parameter  can  be 
used  for  a  wide  range  of  problems. 


Recommendations  for  Three-Dimensional  Modeling 

263.  Unlike  the  two-dimensional  models,  there  are  no  three- 
dimensional  models  that  can  economically  be  applied  for  long-term  reser¬ 
voir  simulations.  This  is  because  all  of  the  models  are  explicit,  and 
thus  excessive  computing  time  is  required.  Imposing  the  rigid-lid 
approximation  removes  the  Courant  condition  on  the  time  step,  but  results 
in  a  Poisson  equation  for  the  pressure  that  must  be  solved,  which  can 
be  costly  in  itself.  Making  the  hydrostatic  pressure  assumption  helps 
in  that  only  a  2-D  Poisson  equation  rather  than  a  full  3-D  Poisson 
equation  must  be  solved.  However,  it  is  not  believed  that  the  rigid-lid 
approximation  is  appropriate  for  models  to  be  used  over  flooding  cycles; 
therefore,  new  solution  techniques  that  allow  for  a  free  surface  but 
remove  the  speed  of  a  gravity  wave  from  the  stability  criteria  must  be 
devised. 


155 


26k.  The  Spraggs  and  Street  (2975)  model  currently  solves  for 
the  free  surface  implicitly,  but  does  not  implicitly  couple  the  internal 
flow  to  these  computations.  It  is  recommended  that  a  coupling  similar 
to  that  in  Edinger  and  Buchak's  (1979)  work,  but  now  in  two  dimensions, 
be  investigated  during  the  next  year.  If  this  can  be  accomplished  in 
an  efficient  manner,  long-term  three-dimensional  free  surface  hydrody¬ 
namic  modeling  will  have  taken  a  giant  step  forward. 
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APPENDIX  A:  NUMERICAL  RESULTS  FROM  APPLICATION  OF 
3-D  PLUME  TO  GRH  FLUME 


1.  The  nonhydro static  version  of  the  Waldrop-Tatom  (1976)*  model 
was  used  to  generate  a  numerical  solution  to  the  density  underflow  problem 
in  the  GRH  flume.  Values  of  various  input  parameters  are  presented  in 
Table  Al.  In  an  effort  to  increase  the  computational  time  step,  and 
thereby  reduce  the  number  of  steps  required,  the  free  surface  was  ini¬ 
tially  assigned  zero  slope  and  was  "frozen"  for  all  subsequent  compu¬ 
tations  . 

2.  With  a  time  step  of  0.5  sec,  the  numerical  solution  was  marched 
forward  for  U000  time  steps  corresponding  to  2000  sec  of  real  time  in 
the  flume.  Outputs  of  velocity,  in  terms  of  the  u  ,  v  ,  and  w  com¬ 
ponents,  and  temperature,  as  taken  from  Tat on  and  Smith  ( 1979a )»  are 
presented  at  0,  1000,  and  2000  sec  in  Tables  A2  through  A13. 


*  References  used  in  the  appendixes  of  this  report  are  listed  in  the  Ref- 
eri  nces  section  at  the  end  of  the  main  text. 


Al 


Table  A1 
3-D  Flume  Input 


I 


Parameter 


•  aluc- 


Inlet  volumetric  flow  rate 
Exit  volumetric  flow  ratu 
Inlet  area 
Outlet  area 
Inlet  velocity 
Outlet  velocity 
Initial  ambient  velocity 
Inlet  temperature 
Initial  ambient  temperature 
Inlet  density 
Initial  ambient  density 
Inlet  equivalent  diameter 
Inlet  kinematic  viscosity 
Inlet  Reynolds  number 
Chezy  coefficient 
Fanning  friction  factor 
Water  depth  at  inlet 
Water  depth  at  outlet 


0.0022';  ft  /sec 
0.0222c  ftf'sec 
0.0  iv' 

0.25  ft' 

(j.01i'45G  ft /sec 
0.08912  ft /sec 
0  ft /sec 
6l.97°;7 
70.79°? 

62.35o  Ibm/ft3 
62.298  lbm/ftJ 

O.67  ft 

1.19  •  10"v  ft2/see 
2196 

13.75  ft"  "/sec 
o.or'.r,' 

1  ft 
3  ft 
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Table  A3 

Region  3  Predicted  Velocity  and  Temperature  Distribution  in  Vertical  Plane  (i  =  2)  at 

Time  =  1000  sec.  Components  u  and  v 
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X-VELOCITY  COMPONENT,  U  (FT/SEC)  Y-VELOCITY  COMPONENT,  V  (FT/SEC) 
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APPENDIX  B:  RESULTS  Or  VECTORIZATION  OF  3-D  PLUME 


1.  A  series  of  benchmark  runs  were  carried  out  to  compare  the 

performance  of  the  scalar  and  vector  versions  of  3-D  PLUME  (Tatom  and 
Smith  1079b) .  Such  runs  were  all  carried  out  with  the  K-Compiler  on 
the  WES  ASC  Computer  and  consisted  of  two  types :  (a)  numerical  con¬ 

sistency  and  (b)  computational  speed. 

Numerical  Consistency  Results 

2.  Benchmark  runs  concerned  with  numerical  consistency  between 
the  scalar  and  vector  programs  were  relati\ ely  short  and  included  print¬ 
outs  of  variables  not  include!  in  the  normal  output  of  either  program . 
Such  printouts  initially  revealed  a  series  of  minor  discrepancies  in  the 
vector  version  of  the  program  and  also  one  previously  undetected  dis¬ 
crepancy  in  the  scalar  version.  After  such  discrepancies  were  corrected, 
numerical  consistency  within  one  percent  was  achieved. 

Computational  Speed  Results 

3.  Two  types  of  timing  runs  were  carried  out  with  both  the  vector 
and  scalar  versions  of  the  program.  The  first  type  consisted  of  only 
the  initialization  computations  plus  one  time  step  computation  and  was 
designed  to  provide  a  measure  of  initialization  time.  The  second  type 
of  run  extended  for  2000  time  steps  and  was  designed  to  provide  a  mea¬ 
sure  of  tiie  computation  time  associated  with  each  time  step.  The  re¬ 
sults  of  these  timing  runs  (including  compilation  time)  are  summarized 

fable  Bl. 

!|.  As  indicated  in  Table  Bl,  the  compilation  time  for  the  vector 
version  of  4  v’e  program  was  approximately  three  times  the  scalar  compila¬ 
tion  time,  as  would  be  expected  because  of  the  additional  optimization 
procedures  involved  in  vector  compilation.  The  initialization  times 
were  essentially  equal.  The  time  required  for  2000  time  steps  with  the 
vector  version  was  approximately  LL  percent  of  the  i-orresponding  time 


Bl 


with  the  scalar  version.  It  is  important  to  note  that  in  carrying  out 
the  timing  comparison  the  amount  of  central  memory  (versus  ext en  led 
memory),  as  shown  in  Table  B1 ,  was  not  the  same  for  all  runs.  For  the 
2000-time  step  vector  run,  no  central  memory  was  used;  while  for  tne 
corresponding  scalar  run  central  memory  comprised  11  percent  of  the 
total  memory.  Because  fetch  times  associated  with  central  memory  are 
smaller  than  the  fetch  time?  with  extended  memory  by  a  factor  of  appro 
imately  6,  the  actual  improvement  in  performance  I  the  vector  version 
over  the  scalar  version  is  sorawnat  greater  than  the-  56-percent  reduc¬ 
tion  in  computat ion  time  n  *ei.  With  an  e  ;jul  portion  of  contra!  mem¬ 
ory,  the  vector  version  should  be  approximately  three  times  as  fast  as 
the  scalar  version.  For  a  production  run  consisting  of  !000  time  step 
tiie  computation  time  would  thus  be  reduced  from  2k 30  sec  to  300  sec. 


1 


Table  B1 

Comparison  of  Computation  Time* 


Scalar 

Vector 

Item 

Central 

Memory 

words 

Extended 

Memory 

words 

Time 

sec 

Central 

Memory 

words 

Extended 

Memory 

words 

Time 

sec 

Compilation 

0 

184,320 

576.63 

0 

184,320 

1,799.89 

Initialization 
and  1  step 

4,096 

106,496 

3.74 

4,096 

131,072 

3.26 

Computation  for 
2000  steps 

12,288 

98,304 

1,200.59 

0 

143,360 

530.14 

f 


*  Tatom  and  Smith  1979b. 
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APPENDIX  C :  NOTATION 


a,b,c 


A 

cv 

*H 


A. 

1 


B 


B. 

1 

B 

o 

C 


C* 


CSHE 


f 

f 


F,f 

F 

s 

g>gi 

G 


Constants  in  expression  for  surface  heat  flux 
Area  of  discrete  element 
Horizontal  eddy  diffusivity 
Cross-sectional  area 

Eddy  diffusivity  tensor  from  time  averaging 
Eddy  diffusivity  tensor  from  spatial  averaging 
Vertical  eddy  diffusivity 
Width 

Surface  width 

Width  of  opening  at  the  dam 

Chezy  coefficient;  constant  in  expression  for  surface  wind 
stress;  phase  velocity 

Constant  in  expression  for  surface  wind  stress 

Sum  of  eddy  diffusivities  due  to  time-  and  spatial-averaging 

Coefficient  of  surface  heat  exchange 

Discrete  element  volume 

Volume  of  a  discrete  element 

Diffusivity  tensor 

Equilibrium  temperature 

Nondimens ional  friction  factor 

Force  vector 

Arbitrary  variables 

Smoothed  solution  in  leapfrog  scheme 

Acceleration  due  to  gravity 

Volumetric  flow  rate 


Cl 


df 

h 

n 

H 


H. 

1 


i  >  J 

J 

k 


L  ,L 
x’  y 


z 

n 


P 

a 


P 

s 


Q 


R 

e 

R. 

1 

R. 

1 

c 


Water  depth 

Height  of  density  underflow 

Rate  of  surface  heat  exchange 

Height  of  opening  at  the  dam;  water  depth 

Surface  elevation 

Unit  vectors 

Boundary  point 

Coefficient  in  expression  for  bottom  friction;  diameter  of 
average  bottom  roughness 

Constant  =  0.10 

Length  scales 

Reference  depth;  length  scale 
Unit  normal  vector  to  the  surface 

Components  of  outward  unit  normal  vector  to  the  surface 
Pressure;  function  to  control  coordinate  spacing 
Time-averaged  pressure 

Time-averaged  and  spatially  averaged  press 

Atmospheric  pressure 

Hydrostatic  pressure 

Reduced  pressure 

Surface  pressure 

Surface  heat  flux 

Function  to  control  coordinate  spacing;  discharge  through 
the  dam 

Reynolds  number 
Richardson  number 
Critical  Richardson  number 


C2 


s 


s.  ,,s 

1  j  ran 
t 

— x 

t 

T 


T" 


T 


T' 


T 


T 

s 


u,v,w 


u. 

1 


u. 

1 


u! 

i 


U 

U 

o 

UWIND 

v 


w 


a 


x,y,z 


x . 

l 


Salinity,  slope  of  reservoir  bottom 

Rate  of  strain  tensor 

Time 

Stress  force  vector 
Temperature 

Deviation  between  instantaneous  and  time-averaged  temperature 
Time-averaged  temperature 

Difference  between  time-averaged  temperature  and  time-  and 
space-averaged  temperature 

Time-averaged  and  spatially  averaged  temperature 

Surface  temperature 

Velocity  components 

Tensor  notation  for  velocity 

Time-averaged  velocity 

Deviation  between  instantaneous  velocity  and  time-averaged 
velocity 

Time-  and  space-averaged  velocity 

Deviation  between  time-averaged  velocity  and  time-  and 
space-averaged  velocity 

Average  velocity  of  density  underflow 

Outlet  velocity 

Velocity  of  the  wind 

Velocity  vector 

Wind  speed 

Cartesian  coordinates 

Tensor  notation  of  spatial  coordinates 
Elevation  of  reservoir  bottom 


C3 


a 


P 

0 


k,£ 

Ap 


At , Ax . 


At 


6.  . 
ij 


e . 


e :  . 
ij 


'ijk 

G 

V 


C,n 


p 

p 

p 


p 


a 


a 


T,T 

o 


n 


Arbitrary  variable;  constant  in  an  expression  for  vertical 
eddy  coefficient  dissipative  coefficient 

Turbulent  Prandtl  number 

Momentum  correction  factor;  constant  in  an  expression  for 
vertical  eddy  coefficient 

Sum  of  ekA  ,  ,  and  p 

Change  in  water  density 

Time  and  spatial  steps 

Time  step  restricted  by  Courant  condition 
Kronecker  delta 
Horizontal  eddy  viscosity 

Eddy  viscosity  tensor  as  a  result  of  time  averaging 

Eddy  viscosity  tensor  as  a  result  of  space  averaging 

Cyclic  tensor 

Vertical  eddy  viscosity 

Water  surface  elevation;  vorticity 

Molecular  eddy  viscosity 

Nonor+.hogonal  curvilinear  coordinates 

Water  density 

Time-averaged  water  density 

Time-averaged  and  spatially  averaged  water  density 
Air  density 

Reference  water  density 
Transformed  vertical  coordinate 
Bottom  shear  stress 
Laminar  stress  tensor 
iiormal  internal  stress  at  the  surface 

Cl 


T 


t 


T  T 
w*  WIND 


Cl) 

V 


Tangential  internal  stress  at  the  surface 
Wind  shear  stress 

Arbitrary  variable;  angle  of  wind  with  reservoir  axis 

Scaling  parameter  in  an  expression  for  the  eddy  viscosity 
tensor 

Coriolis  parameter 
Vorticity 

3/3x  i  +  3/3yj  +  3/3zk 
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In  accordance  with  letter  from  DAHN-RDC,  DAEN-ASl' dated 
22  July  1977,  Subject:  Facsimile  Catalog  Cards  for 
Laboratory  Technical  Publications,  a  facsimile  catalog 
card  in  Library  of  Congress  MARC  format  is  reproduced 
below . 


Johns 'n,  Hilly  H 

A  rev 3  ew  of  nur: er  i  ■  •  a  1  reservoir  hydrodynamic  model  1  nr  /  Ly 
Billy  H.  Jenin':  ei .  (llyiruul  ies  Laboratory.  U .  S .  Arm;, 

Eii.;:nwr  Waterway.-  Experiment  Cation)  ;  pro  i  f  r  Office, 
L’hl-f  >f  Eufiiu-,.rr,  II. P.  Army,  under  r.VQOC  Wore  Unit.  Mo.  !/■,.!.  - 
VicT./Jury ,  Mir;-,  ;  L1 .11.  Arr.y  Knrineor  Waterways  Experiment 
etat.  n  ;  C:  r  i  nr.fi  ■  •  1  i ,  Vu.  :  available  from  STIC,  I’J.'il . 

lr-  »  [C.’l  :.  .  ill.  ?'l  err. .  —  (Technical  report  /  U.C.  Army 
F.nriueer  Waterway;  :ix:  ri merit  Cation  ;  K-oj-I) 
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I.  United  '..Later.  Army.  Corps  of  Engineers.  Office  >f  the 
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